Powers D.L. Boundary Value Problems and Partial Differential Equations

Подождите немного. Документ загружается.

398 Chapter 7 Numerical Methods

x:0.00.20.40.60.81.0

u(x):1.00.643 0.302 −0.026 −0.406 −1.0

Tab le 1 Approximate solution of Eqs. (1) and (2)

Differential equation Boundary condition

u(x) → u

u(0) → u

(x) →

i+1

−2u

i−1

(x)

(0) →

−u

−1

2 x

(x) →

i+1

−u

i−1

2 x

u(1) → u

f (x) → f (x

)

(1) →

n+1

−u

n−1

2 x

Tab le 2 Constructing replacement equations

First, the values of x for the table will be uniformly spaced across the interval

0 ≤x ≤ 1, which we assume to be the interval of the boundary value problem

=i x,x =

These are called meshpoints. The numbers approximating the values of u are

∼

u(x

), i =0, 1,...,n.

These numbers are required to satisfy a set of equations obtained from the

boundary value problem by making the replacements shown in Table 2. The

entry f (x) refers to any coefﬁcient or inhomogeneity in the differential equa-

tion.

Example.

The boundary value problem in Eqs. (1) and (2) would be replaced by the

algebraic equations

i+1

−2u

i−1

(x)

−12x

=−1, i = 1, 2,...,n − 1, (3)

=1, u

=−1. (4)

Equation (3) holds for i = 1,...,n − 1, so the unknowns u

,...,u

n−1

would

be determined by this set of equations. The equations become speciﬁc when we

choose n.Letustaken = 5, so x = 1/5, and the four (i = 1, 2, 3, 4) versions

Chapter 7 Numerical Methods 399

of Eq. (3) are

25(u

−2u

) −

=−1,

25(u

−2u

) −

=−1,

25(u

−2u

) −

=−1,

25(u

−2u

) −

=−1.

(5)

When we use the boundary conditions

=1, u

=−1(6)

and collect coefﬁcients, the foregoing equations become

−52.4u

+ 25u

=−26

25u

− 54.8u

+ 25u

=−1

25u

− 57.2u

+ 25u

=−1

25u

− 59.6u

= 24.

(7)

This system of four simultaneous equations can be solved manually by elimi-

nation or by software. The result will be a set of numbers giving the approxi-

mate values of u at the points x

= 0.2,...,x

= 0.8. The numbers in Table 1

were obtained by a similar process, but using n = 100 instead of n = 5.



Example.

To see how to handle derivative boundary conditions, we solve the problem

−10u = f (x), 0 < x < 1, (8)

u(0) =1,

(1) =−1, (9)

f (x) =







0, 0 < x <

−50, x =

−100,

< x < 1.

The replacement equations for this problem are easily obtained by using Ta-

ble 2. They are

i+1

−2u

i−1

(x)

−10u

=f (x

), (10)

=1,

n+1

−u

n−1

2 x

=−1. (11)

400 Chapter 7 Numerical Methods

We need to know u

, u

,...,u

. The derivative boundary condition at x = 1

forces us to include u

n+1

among the unknowns, so we will need to use Eq. (10)

for i =1, 2,...,n in order to have enough equations to ﬁnd all the unknowns.

Since we have no use for u

n+1

, the usual practice is to solve the boundary-

condition replacement for u

n+1

n−1

−2 x, (12)

and then to use this expression in the version of Eq. (10) that corresponds to

i =n. Thus, the equation

n+1

−2u

n−1

(x)

−10u

=f (x

)

is combined with Eq. (12) to get

n−1

−2 x −2u

(x)

−10u

=f (x

). (13)

Then Eq. (10) for i = 1,...,n−1 and Eq. (13) give n equations that determine

unknowns u

, u

,...,u

To b e s p e c i ﬁ c , l e t u s t a k e n =4, so x =1/4. The three (i = 1, 2, 3) versions

of Eq. (10) are

16(u

−2u

) −10u

=0 (i = 1),

16(u

−2u

) −10u

=−50 (i = 2),

16(u

−2u

) −10u

=−100 (i = 3)

and Eq. (13) adapted to n = 4is



−

−2u



−10u

=−100.

When these equations are cleaned up and the boundary condition u

=1is

applied, the result is the following system of four equations:

−42u

+ 16u

=−16

16u

− 42u

+ 16u

=−50

16u

− 42u

+ 16u

=−100

32u

− 42u

=−92.

(14)

In Table 3 are shown the values of u

obtained by solving Eq. (14) and also

more exact values found by using n = 100.



Elimination is not the only way to get the solution of a system like Eqs. (7)

or (14). An alternative is an iterative method, which generates a sequence of

approximate solutions. For one such method, we solve algebraically the ith

Chapter 7 Numerical Methods 401

x:00.25 0.50.75 1

u (n = 4): 1 2.174 4.707 7.057 7.567

u (n = 100): 1 2.155 4.729 7.125 7.629

Tab le 3 Approximate solution of Eqs. (8) and (9)

equation for the ith unknown. In the resulting set of equations there are “cir-

cular references”: The equation for u

refers to u

and u

, while the equations

for these refer to u

, etc. We may start with some guessed values for the u’s,

feed them through the equations to get improved values for the u’s, and repeat

the process until the values settle down. This method requires a lot of arith-

metic but no strategy, while elimination is just the reverse. It may also work

with nonlinear equations, where elimination cannot.

So far, we have given no justiﬁcation for the procedure of constructing re-

placement equations. The explanation is not difﬁcult; it depends on the fact

that certain difference quotients approximate derivatives. If u(x) is a function

with several derivatives, then

u(x

i+1

) −u(x

i−1

)

2 x

= u



) +

(x)

(3)

(¯x

), (15)

u(x

i+1

) −2u(x

) +u(x

i−1

)

(x)

= u



) +

(x)

(4)

(

¯x

), (16)

where ¯x

and

¯x

are points near x

Now suppose that u(x) is the solution of the boundary value problem

+k(x)

+p(x )u(x) = f (x), 0 < x < 1, (17)

αu(0) −α



(0) = a,βu(1) +β



(1) =b. (18)

If u(x) has enough derivatives, then at any point x

= i x it satisﬁes the dif-

ferential equation (17) and thus also satisﬁes the equation

u(x

i+1

) −2u(x

) +u(x

i−1

)

(x)

+k(x

)

u(x

i+1

) −u(x

i−1

)

2 x

+p(x

)u(x

) = f (x

) +δ

(19)

where

(x)

(4)

(

¯x

) +k(x

)

(x)

(3)

(¯x

Because δ

is proportional to (x)

, it is very small when x is small.

The replacement equation for Eq. (17) is, according to Table 2,

i+1

−2u

i−1

(x)

+k(x

)

i+1

−u

i−1

2 x

+p(x

=f (x

). (20)

402 Chapter 7 Numerical Methods

Thus, the values of u at x

, x

,...,x

, which satisfy Eq. (19) exactly, will nearly

satisfy Eq. (20); vice versa, the numbers u

, u

,...,u

, which satisfy the re-

placement equations (20), nearly satisfy Eq. (19). It can be proved that the

calculated numbers u

, u

,...,u

do indeed approach the appropriate values

of u(x

) as x approaches 0 (under continuity and other conditions on k(x),

p(x), f (x)).

EXERCISES

1. Set up and solve replacement equations with n = 4 for the problem

=−1, 0 < x < 1,

u(0) = 0, u(1) =1.

2. Solve the problem of Exercise 1 analytically. On the basis of Eqs. (15)

and (16), explain why the numerical solution agrees exactly with the ana-

lytical solution.

3. Set up and solve replacement equations with n = 4 for the problem

−u =−2x, 0 < x < 1,

u(0) =0, u(1) = 1.

4. Solve the problem in Exercise 3 analytically, and compare the numerical

results with the true solution.

5. Set up and solve replacement equations with n = 4 for the problem

=x, 0 < x < 1,

u(0) −

(0) = 1, u(1) = 0.

6. Solve the problem in Exercise 5 analytically, and compare the numerical

results with the true solution.

7. Set up and solve replacement equations for the problem

+10u = 0, 0 < x < 1,

u(0) =0, u(1) =−1.

7.2 Heat Problems 403

Use n = 3andn = 4. Sketch the results and explain why they vary so

much.

In Exercises 8–11, set up and solve replacement equations for the problem

stated and the given value of n. If a computer is available, also solve for n twice

as large, and compare results.

−32xu = 0, 0 < x < 1,

u(0) =0, u(1) = 1(n = 4).

−25u =−25, 0 < x < 1,

u(0) =2, u(1) +u



(1) =1(n = 5).

10.

1 +x

=−1, 0 < x < 1,

u(0) =0, u(1) = 0(n = 3).

11.

−u =−x,

(0) =0, u(1) = 1(n = 3).

12. Use the Taylor series expansion

u(x + h) = u(x) + hu



(x) +



(x) +

(3)

(x) +

(4)

(x) +···

with x =x

and h =±x (x

+x =x

i+1

, x

−x =x

i−1

)toobtainrep-

resentations similar to Eqs. (15) and (16).

7.2 Heat Problems

In heat problems, we have two independent variables x and t,assumedtobe

in the range 0 < x < 1, 0 < t. A table for a function u(x, t) should give values

at equally spaced points and times,

=i x, t

=m t,

for i = 0, 1,...,n and m = 0, 1,....Here,x = 1/n, as before. We will use a

subscript to denote position and a number in parentheses to denote the time

level for the approximation to the solution of a problem. That is,

(m)

∼

u(x

, t

404 Chapter 7 Numerical Methods

The spatial derivatives in a heat problem will be replaced by difference quo-

tients, as before:

∂

∂x

, t

) →

i+1

(m) −2u

(m) +u

i−1

(m)

(x)

, (1)

∂u

∂x

, t

) →

i+1

(m) −u

i−1

(m)

2 x

. (2)

For the time derivative, there are several possible replacements. We limit our-

selves to the forward difference

∂u

∂t

, t

) →

(m +1) −u

(m)

t

, (3)

which will yield explicit formulas for computing.

Now, to solve numerically the simple heat problem

∂

∂x

∂u

∂t

, 0 < x < 1, 0 < t, (4)

u(0, t) = 0, u(1, t) = 0, 0 < t, (5)

u(x, 0) = f (x), 0 < x < 1, (6)

we set up replacement equations according to Eqs. (1)–(3). Those equations

are

i−1

(m) −2u

(m) +u

i+1

(m)

(x)

(m +1) −u

(m)

t

, (7)

supposed valid for i = 1, 2,...,n −1andm = 0, 1, 2,....

The point of using a forward difference for the time derivative is that these

equations may be solved for u

(m +1):

(m +1) = ru

i−1

(m) +(1 −2r)u

(m) +ru

i+1

(m), (8)

where r = t/(x)

.Thuseachu

(m +1) is calculated from u’s at the preced-

ing time level. Because the initial condition gives each u

(0),thevaluesofthe

u’s at time level 1 can be calculated by setting m =0 in Eq. (8):

(1) =ru

i−1

(0) +(1 −2r)u

(0) +ru

i+1

(0).

Then the values of the u’s at time level 2 can be found from these, and so on

into the future. Of course, r has to be given a numerical value ﬁrst, by choosing

x and t.

It is convenient to display the numerical values of u

(m) in a table, making

columns correspond to different meshpoints x

, x

,...,x

and making rows

correspond to the different time levels t

, t

,....SeeTable4.

7.2 Heat Problems 405

m 01234

0 00.25 0.50.75 1

1 0 0.25 0.50.75 0

2 0 0.25 0.50.25 0

3 0 0.25 0.25 0.25 0

4 0 0.125 0.25 0.125 0

5 0 0.125 0.125 0.125 0

Tab le 4 Numerical solution of Eqs. (4)–(6)

Example.

Solve Eqs. (4)–(6) with x = 1/4andr = 1/2, making t =1/32. The equa-

tions giving the u’s at time level m +1are

(m +1) =



(m) +u

(m)



(m +1) =



(m) +u

(m)



(m +1) =



(m) +u

(m)



(9)

Recall that the boundary conditions of this problem specify u

(m) = 0and

(m) = 0 for m = 1, 2, 3,.... Thus we ﬁll in the columns of the table that

correspond to points x

and x

with 0’s (shown in italics in Table 4). Also

the initial condition speciﬁes u

(0) = f (x

),sothetoprowofthetablecanbe

ﬁlled. In this example we take f (x) = x, and the corresponding values appear

in italics in the top row of Table 4.

The initial condition, u(x, 0) = x,0< x < 1, suggests that u(1, 0) should

be 1, while the boundary condition suggests that it should be 0. In fact, nei-

ther condition speciﬁes u(1, 0), nor is there a hard and fast rule telling what to

do in case of conﬂict. Fortunately, it does not matter much, either. (See Exer-

cise 1.)



Stability

Thechoicewemadeofr = 1/2 in the Example seems natural, perhaps, because

it simpliﬁes the computation. It might also seem desirable to take a larger value

of r (signifying a larger time step) to get into the future more rapidly. For

example, with r = 1(t = 1/16) the replacement equations take the form

(m +1) = u

i−1

(m) −u

(m) +u

i+1

(m).

In Table 5 are values of u

(m) computed from this formula. No one can believe

that these wildly ﬂuctuating values approximate the solution to the heat prob-

406 Chapter 7 Numerical Methods

m 01 2 34

0 00.25 0.50 0.75 1

1 0 0.25 0.50 0.75 0

2 0 0.25 0.50 −0.25 0

3 0 0.25 −0.50 0.75 0

4 0 −0.75 1.50 −1.25 0

5 0 2.25 −3.50 2.75 0

Tab le 5 Unstable solution

lem in any sense. Indeed, they suffer from numerical instability due to using a

time step too long relative to the mesh size. The analysis of instability requires

familiarity with matrix theory, but there are some simple rules of thumb that

guarantee stability.

First, write out the equations for each u

(m +1):

(m +1) = a

i−1

(m) +b

(m) +c

i+1

(m).

The coefﬁcients must satisfy two conditions

1. No coefﬁcient may be negative.

2. The sum of the coefﬁcients is not greater than 1.

In the example, the replacement equations were

(m +1) = ru

(m) +(1 −2r)u

(m) +ru

(m),

(m +1) = ru

(m) +(1 −2r)u

(m) +ru

(m),

(m +1) = ru

(m) +(1 −2r)u

(m) +ru

(m).

The second requirement is satisﬁed automatically, because r +(1−2r) +r = 1.

But the ﬁrst condition is satisﬁed only for r ≤ 1/2. Thus the ﬁrst choice of

r =1/2 corresponded to the longest stable time step.

Example.

Different problems give different maximum values for r. For the heat conduc-

tion problem

∂

∂x

∂u

∂t

, 0 < x < 1, 0 < t, (10)

u(0, t) = 1,

∂u

∂x

(1, t) +γ u(1, t) = 0, 0 < t, (11)

u(x, 0) = 0, 0 < x < 1 (12)

7.2 Heat Problems 407

the replacement equations are found to be (for n = 4)

(m +1) = ru

(m) +(1 −2r)u

(m) +ru

(m),

(m +1) = ru

(m) +(1 −2r)u

(m) +ru

(m),

(m +1) = ru

(m) +(1 −2r)u

(m) +ru

(m),

(m +1) = 2ru

(m) +



1 −2r −

rγ



(m).

(13)

(Remember that u(1, t), corresponding to u

, is an unknown. The boundary

condition has been incorporated into the equation for u

(m + 1).) Again, the

second stability requirement is satisﬁed automatically; but the ﬁrst rule re-

quires that

1 −2r −

rγ ≥ 0orr ≤

2 +

. (14)



EXERCISES

1. Solve Eqs. (4)–(6) numerically with f (x) = x,asinthetext(x = 1/4,

r =1/2), but take u

(0) = 0. Compare your results with Table 4.

2. Solve Eqs. (4)–(6) numerically with f (x) = x, x = 1/4, u

(0) = 1, as in

the text, but use r = 1/4. Compare your results with Table 4. Be sure to

compare results at corresponding times.

3. For the problem in Eqs. (10)–(12), ﬁnd the longest stable time step when

γ = 1, and compute the numerical solution with the corresponding value

of r.

4. Solve the problem in Eqs. (10)–(12) with x = 1/4, r = 1/2andγ = 0, for

m up to 5.

For each problem in the following exercises, set up the replacement equations

for n = 4, compute the longest stable time step, and calculate the numerical

solution for a few values of m.

∂

∂x

∂u

∂t

, u(0, t) = u(1, t) = t, u(x, 0) =0.

∂

∂x

−u =

∂u

∂t

, u(0, t) = u(1, t) = 1, u(x, 0) = 0.

∂

∂x

∂u

∂t

−1, u(0, t) = u(1, t) = 0, u(x, 0) = 0.