Chapter 4 Elementary Examples of Difference Schemes

Chapter 4 Elementary Examples of Difference Schemes

71 Part 2 DIFFERBNCE SCHEMES FOR ORDINARY DIFFERENTIAL EQUATIONS Part 2 of this book is devoted to the construction and the study of difference sche...

421KB Sizes 2 Downloads 118 Views

71

Part 2 DIFFERBNCE SCHEMES FOR ORDINARY DIFFERENTIAL EQUATIONS

Part 2 of this book is devoted to the construction and the study of difference schemes for ordinary differential equations. In the course of this study we introduce the concepts of convergence, approximation and stability, basic in the theory of difference schemes and general in character. Familiarity with these concepts, acquired in connection with ordinary differential equations, will permit us later, in the study of difference schemes for partial differential equations, to concentrate on the numerous peculiarities and difficulties characteristic of this most variegated class of problems.

Chapter 4 Elementary Examples of Difference Schemes In this chapter we consider introductory examples of difference schemes, intended only to give the reader a preliminary acquaintance with the basic concepts of the theory. SS.

The concept of order of accuracy and of approximation

1 . Order of accuracy of a difference scheme. This section is devoted to the question of the convergence of solutions of difference equations, with refinement of the net, to the solutions of the differential equations which they approximate. We limit ourselves, here, to the study of two difference schemes for the numerical solution of the problem

du -+au=O, dx

O(x(1,

u ( 0 ) = b.

Let us begin with the simplest difference scheme, based on the use of the difference equation

We now subdivide the interval [0,11 into steps of length h. It is convenient to take h = 1/N, where N is an integer. The points of

72

Chapter 4

Elementary Examples of Difference Schemes

subdivision will be numbered from left to N. The value of u obtained, via the will be denoted as un. We fix an initial From difference equation (2) one gets the

...,

= (1

u

-

Ah)u

right, so that x, = nh, n = 0 , 1, difference scheme, at point xn value, uo. Suppose that uo = b. relation n-1'

from which we find the solution of Eq. (2) subject to the initial condition uo = b: u

= (1

-

b = (1

Ah)n

-

Ah)

(xn/h)

The ezuct solution of problem (1) has the form u(x) takes on, at point xn, the value u(xn)

=

be-Axn

(3)

b.

= b exp

.

=

[(l

-

Ah)

It

(4)

We now estimate the error in the approximate solution ( 3 ) . this error is S(xn)

(-Ax).

bnh)- e-hn]

At point xn

b.

(5)

We are interested in the rate at which 6(x ) decreases as the number o f subdivision points increases or, equivalently, as one decreases the stepwidth, h ( ; , f [ T , of the difference net. To bring this out we represent

(1

- Ah)

in the form X

X

n i;n- ln(1 (1-fijh = e

-

X

f

Ah)

[-Ah

= e

A2 h2 + -y-+

O(h3)]

+ = e

A2x

-Ax

+ h+e

O(h2)

" + O(h2).

Thus Eq. (3) takes the form -Ax

u

=be

A2xn " + h b r e-Axn

+

O(h2),

so that

6(xn) = hb

A2xn

-Axn

7 e

+

O(h2)

=

O(h),

1

[l

+

O(h2)]

=

$8

73

Order of Accuracy and Approximation

i.e. the error ( 5 ) tends to zero as h + 0, and the magnitude of the error is of the order of the first power of the step-size. On this basis one says that the difference scheme has first-order accuracy (which is not to be confused with the order of the difference equation, defined in 81). Let us now solve problem ( 1 ) with the aid of the difference equation

This is not as simple as it may seem to be at first glance. The problem is that the above scheme is a difference equation of second-order, i.e. it requires the assignment of two initial conditions (u(xo) = 0 and u(x,) = u(h)); while the equation to be integrated, Eq. ( I ) , is an equation of first order, and for it we need only the condition u(0) = b. It is natural, also in the difference scheme, to set uo = b. It isn't clear, however, how one should choose ul. To shed some light on this question we use the explicit form of the solution of Eq. ( 7 ) (see 83 Eq. ( 6 ) ) :

where 91

f i - Ah

=

q2 = (-1)(1

+

=

Ah

1

-

Ah

A2h2 2

-I--+

O(h4), (9 1

+1 A2h2 ) + O(h4).

The Taylor expansions, (91, of the roots of the characteristic equation We carry allow one to develop an approximate representation of qn and 9;. out a detailed derivation of such a representation for

in: 1

Since Ln(1

+

z) = z

1

Therefore

-

z2/2

+

z3/3

A2h2 Ah + T +

+ O(z4),

O(h4)]=

-Ah

+-A3h3 6 + O(h4).

74

Chapter 4

Elementary Examples of Difference Schemes X

q;1 = e

[-Ah

A3 h3 +7 +

O(h4)]

[

-AX, = e

1

+ h2

+] A3x

+

O(h3).

(10)

We will not carry out the completely analogous computation for q;, go directly to the result: qn2 = (-l)”ekun

+

O(h2)

hut

(11)

Putting the approximate expressions for q y and q; into E q . ( 8 ) we get

+

q l u 0 - u1 ( - l ) n [ekn 42 - 91

O(h2)].

All further conclusions will be obtained through study of this expression. We note that if the coefficient (q2uo - u , ) / ( q 2 - q ) tends to the 1, q l ) , on finite limit b as h + 0 then the first term, (q2u0 - u l ) q l / ( q 2 the right-hand side of E q . ( 1 2 ) tends to the desired solution of problem

-

(I)., Since

i . e . does’not converge to a dePinite limit, then to guarantee convergence to a limit, as h 0 , of the second tern on the right-hand side of (12), -f

+

‘lUo- u1 (-l)n[ehn 42 - 41

O(h2)],

it is necessary to require that the expression (qluo to zero as h + 0. Let us, then summarize what has been said. So that the solution of the difference equation U(X

+

h)

- U(X

2h

-

h)

+

-

u,)/(q2

-

ql)

Au(x) = 0

should converge to the solution u = b exp (-Ax) of the boundary-value problem (l), it is necessary that the conditions

tend

75

Order of Accuracy and Approximation 4lU0

-

u1

-

-

42u0

0,

u1

+

b.

ql 42 - 41 be satisfied. Recall, further, that we chose to set uo equal to b. Condition ( 1 4 ) gives us a hint as to how we can assign u l . It turns out that it is sufficient that u1 + uo as h + 0 . In fact q 1 + + l , and q2 + -1 a s h + 0 and, therefore, as h + 0 42

q1u0 42

+

- 1 - 91

-

q2uo +

0,

-

42

+

91

b.

2. Speed of convergence of the solution of the difference equation. We now go on to a study of the speed of convergence for different specific choices of u ufh). 1 To determine u(h) it is natural to make use of the Taylor series expansion of the solution of the differential equation u' + Au = 0 . Using the fact that u' = -Au, we rewrite the Taylor series expression thus: J

u(x,)

=

u(0)

-

hAu(0)

+

-

O(h2) = u(O)(1

+

Ah)

O(h2).

This equation is satisfied by the exact solution of the differential equation. In the approximate solution, limiting ourselves to two terms of this expansion, we can set u

1

= u

(1 0

-

Ah).

If we have decided to take only one term we let

In the first case we commit, in the initial value u l , an error of order h2, in the second -- an error of order h. Let use examine the speed of convergence in each of these two cases, for each of these assignments of initial values. Assume uo =

b,

-

u1 = (1

Ah)b.

(15)

Then (see Eq. ( 9 ) ) 'lUO 42

q2uo

42

-

-

1 91

-

-

[-I

u1

-

[l

- Ah

41

- Ah

+

O(h2)]b

-2

- 'A

h2 -2

+

-

(1

-

Ah)b

=

O(h2)

+ O(h4)]b + ~(h')

(1

-

Ah)b =

b

0(h21

+

I

O(h2).

(16)

76

Chapter 4

Elementary Examples of D i f f e r e n c e Schemes

Returning t o Eq. (12), we e a s i l y come t o t h e c o n c l u s i o n which h a s been o u r goal un = be-AXn

+ O( h2) .

T h i s c o n c l u s i o n may be s t a t e d as f o l l o w s . If t h e i n i t i a l v a l u e u1 is g i v e n c o r r e c t l y t o o r d e r h 2 , t h e n t h e e r r o r i n t h e s o l u t i o n w i l l be o r d e r h 2 , i . e . t h e d i f f e r e n c e scheme w i l l be a c c u r a t e t o second o r d e r .

It can be shown t h a t , even i f we t a k e , f o r u1, i t s e x a c t v a l u e exp (-Axl), a c c u r a c y g r e a t e r t h a n o r d e r h2 c a n n o t be a t t a i n e d i n t h e b s o l u t i o n . We a d v i s e t h e r e a d e r t o prove t h i s a s s e r t i o n as an e x e r c i s e . It

is easy t o show a l s o t h a t i f , f o r UO, we t a k e n o t p r e c i s e l y b , b u t any q u a n t i t y of t h e form b O ( h 2 ) , t h e speed of convergence w i l l s t i l l be second o r d e r . We now proceed t o c o n s i d e r t h e second f o r m u l a t i o n of i n i t i a l c o n d i t i o n s we have set o u t t o s t u d y . Suppose

+

u1 = uo = b. Now

qluO q2

- ul - 91

[l

-

Ah

+ -2

O(h2)lb

+ O(h2)

-

b

=1 2 Ahb

+

o(h2),

and, c o n s e q u e n t l y , q2uo u

n

=

92 = [b

-

-u -

-Ax 1 [e

91

++

(-1)"[i

Ahb

+

Ahb

+ O(h2)] -Ax O(h2)][e

q1u0 92

-

41

(-l)"[ehn

+ O(h2)]

=

+ O(h2)] -

+ O(h2)][ehn +

O(h2)] =

- (-1p eAXn 2

-u

h

+ O(h2).

Thus i f t h e e r r o r i n i n i t i a l v a l u e s is of o r d e r h , t h e n t h e e r r o r i n t h e s o l u t i o n w i l l a l s o be of o r d e r h.

77

Order of Accuracy and Approximation

50

Let us now summarize what has been said. difference scheme examined earlier,

+

U(X

-

h)

-

h)

u(x)

+

U(X

2h

We have seen that the

+ Au(x)

=

0,

as compared to the scheme

+ h) -

U(X

h

= 0,

Au(x)

can give faster convergence, and more precisely convergence with remainder terms of order h2, rather than order h as in the second of these schemes. In order to attain second order accuracy one must, having taken an exact uo, choose a u differing from the exact solution of the differential 1 equation at point x = xo + h by a quantity of order h2. It can be shown that uo also need not be given exactly, but also may contain an error of order h2. The speed of convergence i s not thereby diminished. Refining the initial values up to order h3 and higher does not result in an increase in the accuracy of the solution. If the initial values are given with errors of order h, then the solution will contain an error of this same order. 3. Order of approximation. It is interesting to consider just what it is that renders the scheme

+ h) -

U(X

U(X)

h

+

= 0

Au(x)

less accurate than the scheme U(X

+

h)

-

U(X

2h

- h) + Au(x)

= 0.

These schemes differ in the approximate expressions u(x

+ h) h

u(x)

and

U(X

+ h) 2h

U(X

- h)

used for the derivative, du/dx, at point X. It is natural to assume that in the first scheme the derivative has been replaced by a less accurate expression than in the second. And this is, in fact, true. Let us substitute, for u(x + h) and u(x - h), their Taylor series expansions u(x

+

h)

u(x)

+

u’(x)h

+

u(x

-

h) = u(x)

-

u’(x)h

+ u”(x)

=

Using these expressions we get

u”(x)

h2 F + u”’(x)

h3 r +O(h4],

h2

h3

-

u”’(x)

+

O(h4).

Chapter 4

Elementary Examples of Difference Schemes

78

u(x

+

h, h

u(x

+ h)

-

u(x)

- U(X

-

2h

+ u”(~)

= u’(x)

h)

h

+ u”-(x)

= u’(x)

+ h2

O(h2).

+ 0(h4),

i.e. in the first case we have an approximation to the derivative of only first-order accuracy, and in the second -- of second order. The examples we have considered might lead one to think that the order of convergence of solutions of difference equations can be taken to be equal to the order of approximation of the derivatives in the differential equation. It turns out, however, that in such a very general form, this hypothesis is untrue. On those difference schemes for which it’s validity will be proven it will be necessary to impose an essential restriction -- the requirement of stability. The necessity of this requirement will become clear as we consider the examples in the following section.

s

9. Unstable difference 6chemes 1. Techniques for approximating the derivative. We now again consider difference schemes for the approximate integration of the simplest differential equation u‘ + Au = 0. A s we have already seen, to construct a difference scheme approximating this equation it suffices to replace the derivative, u’, by some sort of approximating difference expression. Thus we have examined schemes in which the derivative u’ was replaced by u(x + h) h

- u(x>

+ h)

u(x

or

- U(X

2h

-

h)

It is clear than any expression of the form !J

u(x

+ h) -

U(X

2h

-

h)

+

(1

- v)

U(X

+ h) h

u(x)

will also approximate u’(x). In fact let us substitute into this expression the Taylor series for u(x + h) and u(x h):

-

u(x

+ h)

U(X

-

-

h)

= U(X)

h) =

U(X)

+

u’(x)h

+

O(h2)

-

u’(x)h

+

O(h2).

,

We then get

v

u(x

+ h) =!J

- u(x

2h [u(x)

+

-

,,) u(x

+ u’(x)h + O(h2)] -

2h

+ h) h

[u(x)

U(X)

- u’(x)h

+

O(h2)]

+

$9

79

Unstable Difference Schemes

Using this sort of approximation for the derivative one can derive a whole family of difference schemes depending on a numerical parameter. These schemes will have the form

To each value of the parameter p corresponds one such scheme. It was the study of those particular schemes for which P = 0 and v = 1 to which $8 was devoted. 2. Example of an unstable difference scheme. We now consider one more scheme of this form, obtained from (1) with v = 4 :

4

U(X

+

h)

-

U(X

-

h)

2h

- 3

+

U(X

h) h

-

U(X)

+ Au(x)

= 0.

(2)

This scheme may be rewritten thus: -2u(x

- h) +

(3

+

-

Ah)u(x)

U(X

+ h)

= 0.

As in the examples considered earlier, we compute the solution interval [0,1], subdivided by the points of the difference net steps, each of length h = 1/N. The coordinate, x of a point n’ is defined a s x = nh = n/N. The solution of the difference equation may be written in form

where q1 and q

2

(2’ ) on the into N equal of the net the explicit

are roots of the characteristic equation -2

+

(3

+ Ah)

q

-

+

2A2h2

q2 = 0 .

Let us compute q1 and q2: + Ah

q1 =

42

=

+

Ah

+

+

6Ah +

+

6Ah

2

2

=

1

-

Ah

+

O(h3), (4)

+

= 2 ( 1 + Ah) + o(h2).

We will use, for qy and qi, the approximate expressions

80

91 q;

Chapter 4

Elementary Examples of Difference Schemes

= [l

-

= [2(1

Ah

+ 0(h2)]"

+ Ah) +

= [l

0(h2)]"

= [2(1

=

+ Ah) + O(h2)]

(Xn/h) = 2

.1

-Ax

(xn/h)

- Ah + O(h2)]

Axn

e

+ O(h),

(xn/h)

(5)

+ O(h)l.

Substituting Eq. (5) into ( 3 ) we get

Before considering what limit u will tend to as h + 0 we must indicate how n we will fix the initial values, uo and ul, of the difference solution. Just as in 58 we will look for a solution satisfying the condition u(0) = b, 2nd take as difference starting values uo = b and u = b(1 - Ah). 1 We substitute these starting values into Eq. ( 6 ) and simplify each term separately. The first and second terms, respectively, take the forms 92uo 92

- u1

-

-Axn [e

41

-

+

[2 [2

+

[I

-

O(h)]b O(h)]

+ O(h)]

=

- (1 - Ah)b - [ 1 - Wh)]

-Ax

Ah + 2A2h2 + O(h3)lb [l + O(h)] - [2 + O(h)]

-Ax

+ O(h)]

= be

- b(1 - Ah)

teAxn

[e

=

-

+

+

O(h),

O(h)]

(Xn/h)

+

2

2A2h2b[eAxn

O(h)]

-

(xn/h)

Thus we get

For x = x = const, as h * 0 the first term of this expression tends n to b exp(-Ax), i.e. to the desired solution. Therefore if the whole expression for un is t o converge to this solution it i s necessary that the second term should go to zero: but as h + 0 this term tends, not to zero, but to infinity. In fact -2A2b exp ( A x ) + O(h) tends to the finite,

59

81

Unstable Difference Schemes

( xn/h) nonvanishing, limit -2A2b exp(Ax) , and h2 2 tends to infinity faster than any positive power of lfh. We have shown that a difference scheme approximating the differential equation can have a solution not converging, as h + 0 , to the solution-of the differential equation. One might think that the fault lies, here, in an insufficiently accurate choice of ul. But we will now show that there will be no convergence even if we take u1 to be exactly equal to the solution of the differential equation at x1 = xo + h, that is if we set u1 = uo exp(-Ah) = b exp (-Ah). Let us begin by simplifying the expressions occurring in Eq. ( 6 ) : qzu0

-

u1 -

92

- 41

91

-

[Z

[2 + O(h) b O(h)J -][1

+

- be-Ah + O(h)]

=

+

o(h)9

92

Substituting these expressions into ( 6 ) we get un

=

[be-Axn

+

O(h)]

-

[S

A2behn

+

O(h)]h2

2(Xn/h)

(7)

The second term on the right-hand side of this equation again tends to infinity, while the first remains bounded. Therefore the whole solution of the difference equation also tends to infinity. The reason that difference scheme (2) doesn't converge as h + 0, as we have seen, is the fact that it can have solutions which grow quickly as the step-size h decreases, even if the starting values are completely reasonable. Such difference schemes are called "unstable". Naturally, they are unsuitable for the numerical solution of differential equations.