DIFFERENCE EQUATIONS

II

DIFFERENCE EQUATIONS

In this chapter we shall review the elementary properties of finite difference equations. These equations can be readily solved by the use of high-speed digital computers. T h i s type of equation arises in practice usually as an approximation to a differential equation whose solution cannot easily be found analytically. T h e prominence of differential equations in science and engineering leads to concentrating our numerical studies on the difference equation.

W e shall be concerned primarily with techniques suited to digital computers, even though they certainly are not a panacea for all problems.

Frequently the exact solution of a problem requires an infinite number of steps. Although a machine can perform an extremely large number of operations, this number is finite. For this reason, the numerical calculation must be terminated at some stage. Usually the termination of a calculation introduces an error, called the truncation error.

Other types of errors arise in the use of computers. A digital computer carries only a finite number of digits, and any number requiring more digits will be represented only approximately. T h e error introduced by carrying a limited number of digits is called the round-off error.

A large portion of the numerical analysts' effort is involved in esti- mating the errors of all types present in a given numerical procedure.

In many cases reliable estimates of the error are unobtainable. Experi- ments or other measures are then necessary to test the accuracy of the approximation. Whenever practical we shall discuss errors; however, our main concern will be with the formulation of methods of solution.

2.1 A Simple Example

A s an introduction to numerical techniques, w e consider the problem of finding the area A under a curve as shown in Fig. 2.1.1. W e assume 53

the curve y(x) is specified as the solution of the following first order differential equation.

( 2 . 1 . 1 ) where f(x) and y(a) are known. T h e analytic solution to this problem is simply

A = j^b dx [J/(*') dx'] + c(b - a), (2.1.2) where the constant c is determined from the condition at χ = a.

T h e elementary definition of a derivative, viz., dy = mU y(x + Δχ)

-

^y(x)

dx AX-*Q Δχ (2.1.3)

leads to a suitable numerical procedure. T h e limit is approximated by the so-called first divided difference:

y(x

+

Δχ) — y(x)

Δχ (2.1.4)

Λ

y(x) _A

a χ x + Ax b

F I G . 2.1.1. A n arbitrary curve the area under which is to be calculated.

T h i s approximation corresponds to replacing the actual curve between χ and χ + Δχ by the secant line. L e t us divide the entire interval from a to b into / equal subintervals of width Δχ. T h e abscissas are numbered from 0 to / starting at χ = a. For the^th abscissa,

Xj = a+j^x, y = 0 , l , . . . , / . (2.1.5) T h e corresponding ordinates are

y(*i) = yt = y(* + βχ). (2.1.6)

2.2 D I F F E R E N C E A N D S U M M A T I O N O P E R A T O R S 55

By Eqs. (2.1.4) and (2.1.6), w e find¹ that

Λ + ι = Λ + / ( * * ) (2.1.7) and

Vi = JO + £ / ( * * ) J = 1 . 2 , / . (2.1.8)

k=0

T o find the approximate area under the curve, w e must estimate the behavior of y(x) at values of χ intermediate between the xi. For simplicity we approximate y(x) by a series of rectangles

y(*) = Vi» * i < * < . (2.1.9)

T h e total area under y(x) in the interval a < # < b is then

^ = !*> + Σ

^{[ λ +}

S/fo) Η1 ^.

(2-1 · 10)

j 2

Λ = y0(b - a) + V X/foXJ*)¹, (2.1.11) an approximate solution whose form is similar to the exact solution

(2.1.2).

T h e present example is typical: divided differences replace derivatives, finite sums replace integrals. T h e independent variable is divided into a finite number of values.

2.2 Difference and Summation Operators

I n order to develop approximations to differential and integral equations, it is convenient to define various difference and summation operators. I n this section w e shall introduce several operators of use and derive a number of relationships between the operators. Further, w e shall show certain properties of the operators which are analogous to properties of differential and integral operators.

L e t h = Axf the spacing between abscissas. T h e first forward difference A is defined by

Ay{x) = y{x + h) - y(x\ (2.2.1)

1 Henceforth the equals sign is used in relations such as (2.1.7), even though the divided difference is only an approximation.

Ac Acf(x)

*U(x) + *(*)]

*if(x)g(x):

[six).

J " [ J »/(*)]

T h e shift operator Ε is defined by

Ey(x) = y(x + A), (2.2.7) the first backward difference V by

Vy(x) = y(x) — y{x — h\ (2.2.2)

and the central difference operator δ by

Sy(x) = y(x + A/2) - y(x - A/2). (2.2.3) T h e first forward difference operator was used in the previous section.

Repeated application of the operators leads to the expression

A"y(x) = J ) ( - ) * k Kn[n ky y[x + ( « - k)h], (2.2.4)

Vy(x) = (- )* Μ ( ηη ]_ ky y(x - kh), (2.2.5) and

8²".rW = X ( - ) * ^fe!(2(f- j ^ - M * + ( » - *)^A1> (²·²·⁶) as may be seen by induction. It has been assumed that the spacing h is constant.

These difference operators have many properties in common with the differentiation operator, as can be seen from T a b l e 2.2.1. T h e relations in the table are easily verified.

T A B L E 2.2.1

PROPERTIES OF THE F O R W A R D DIFFERENCE O P E R A T O R

0 (c = constant) cAftx)

*f(x) + Mx)

f(x)Ag(x) + g(x + h)Af(x) g(x)Af(x) -f(x)Ag(x)

g(x + h)g{x) J »»+»/(*)

2.2 D I F F E R E N C E A N D S U M M A T I O N O P E R A T O R S 57

and the inverse shift operator E'¹ by

E~^ly(x) = y(x - h). (2.2.8)

T h e indefinite summation operator Σ (sometimes denoted Δ~^χ) is defined as follows: if

z{x) = Ay(x), (2.2.9) then

Xz(x)=y(x). (2.2.10) T h e addition of an arbitrary constant is omitted for our purposes.

T h e various operators have many relations among themselves. For examples,

Δ =E-\y (2.2.11)

V = 1 -E~\ (2,2.12)

δ² = J V , (2.2.13)

X J = 1. (2.2.14)

T a b l e 2.2.2 lists various summation formulas.

T A B L E 2.2.2 S U M M A T I O N F O R M U L A S

Σ φ ) = εΣ/(χ) Σ [/(*) + * ( * ) ] = Σ fix) + Σ * ( * )

Σ [f(x)Mx)] = /(*)*(*) -Σg(χ + h)Af(x) Σ [f{x + h)Ag(x)] = f(x)g(X) - ΣΜ(Χ)Δ/(Χ)

T h e factorial series are defined by

= x[x - h][x - 2h]... [χ -(η - 1)Λ], (2.2.15)

= x[x + h][x + 2h] ... [x + {n- \)h\. (2.2.16)

L e t the arguments be denoted by subscripts, viz.,

Xj = x+jh, (2.2.17)

x_j = Λ (2.2.18)

which is analogous to the usual rule for differentiation. A rule analogous to the usual one for integration is found by applying the operator Σ to Eq. (2.2.23).

2*

^{, n ,}

= -FFiyr- <

²

·

²

·

²⁴

>

A series of numbers summed between definite limits is called a definite sum. L e t z(x) be such that

z(x) = Ay(x).

W e observe that

X * ( * ; ) = 3<*i) ³²⁺¹ · (2.2.25) T h e sum formula (2.2.25) is very useful in summing the terms of a

series in n^j. For example, consider the sum of the first η terms of the series in n^z. Since w⁰⁾ is a polynomial of degree j in n> any polynomial of degree j can be expressed in terms of the factorial series up to and includ­

ing degree j. I n the example

n^a) = w,

w(2) = W(w _ ! ) = W2 _ w> (2.2.26)

n<3> = n(n _ i ) (w _ 2) = n^z - 3n² + In,

so that

η» = nw + 3n<²> + w^{( 1 )}. (2.2.27) (2.2.19)

(2.2.20)

(2.2.21)

(2.2.22)

(2.2.23) Application of the operator Δ to xⁱⁿ⁾ shows that

so that

N o t e that

and define

2.3 D I F F E R E N C E E Q U A T I O N S A N D T R U N C A T I O N E R R O R 59

Therefore,

»(4) ~(2)-,ΛΓ+1

?^{w3 =} [ⁱr- ^{+ w ( 3, +} V ] 1 ( 2·²·^{2 8)}

= ^ ^ +(Λ Γ + 1 Γ +( ^ + ^ ( 2.2.29 )

1 )2 • (2.2.30)

2.3 Formation of Difference Equations and Truncation Error

T h e construction of a difference equation from a differential equation is not a unique process. M a n y different difference approximations are possible for a given differential equation. T h e selection of a particular difference relation is usually determined by the nature of the truncation error associated with the approximation.

A s an introductory example, let us develop expressions for the second derivative in terms of the forward, backward, and central difference operators. A uniform spacing h is postulated.

W e assume the function y(x) may be expanded in a Taylor series in the closed interval χ — 2h < χ < χ + 2h. W e have

X* ± k) = * , ) ± hyV) + + .... (2.3.1)

y(* ± 2k) = y(x) ± ( 2 % ' ( * ) + » / ' ( * ) + ... + + .·., (2.3.2) where a prime denotes differentiation. W e have then

y {x + 2h)-2y(X + h)+y(x) ^ y i) +{ x A y) +> ^ + ^ ( 2 < 3 > 3)

or

where the notation 0(h) means the first term neglected is of order h.

Similarly we have

- ^ = ^ + 0(/,). (2.3.5)

However,

8²y(x) _ d²y(x)

h² + 0{h% (2.3.6)

T h e forward and backward differences are accurate to order h\ the central difference, to order h², a result that is intuitively and geometrically obvious. T h e second forward difference approximates the second derivative by taking the difference in slope of the secant lines de and cd in Fig. 2.3.1. T h e backward difference uses the secants ab and be, whereas the central difference uses the lines be and cd. Since the latter pair are centered with respect to the point of interest, we should expect the central difference to be more accurate, as indeed it is.

T o consider higher order approximations and other derivatives, we return to the Taylor series. W e then have

M* + h) = y(X) + h + - » φ"y(X) + . . . + £ ( | f y(x) + (2.3.7) hⁿ /d\ⁿ

or

dX) dX)

y(x). (2.3.8)

x — 2h x—h x+h

F I G . 2.3.1. Approximation of a curve by secant lines.

T h e series in brackets is the expansion for the exponential and hence we have (formally)

d)

E = ^ (^hd i ) ' (2.3.9)

2.3 D I F F E R E N C E E Q U A T I O N S A N D T R U N C A T I O N E R R O R 61 Treating E q . (2.3.9) as an identity, w e may derive expressions for any order derivative in terms of the various difference operators. For example,

- τ = ( i^{l n} )^{Ey {) =x} I^{ln Δ +{} X ) y { x )- ^{( 2}·³·^{1 0)}

Using the logarithmic expansion, w e have dv(x) \ \Λ Δ² J³ ι

- ί ^ = α Η - τ + τ - . ρ ^μ· ^{( 2}·³·^{Π )} For higher order derivatives, w e have

^ = 1L[ l n ( J + l ) ] M * ) , (2-3.12)

=

-L[δ'-\δ«» + J 4 * L+ i L +

...]

y{x

y

( 2. 3 .1 3)

In terms of the backward difference operator, w e have

Μ Χ) = - J l n ( l - V M * ) . (2-3.14) dx h

or

- * 3 ? - ί [ ' + τ + τ + - ] * * <"·'«

T h e higher derivatives are

" 2 ?

^{- Τ}

[" + + ^ f ^ v " +

···] ( " - Ι Ο For the central difference operator, one is usually interested in even powers of the operator. By performing the necessary expansions (see problem 4), w e find

d²y(x) dx and

d²ⁿy(x) _ 1

Ρ ⁼ - Μ ^{δ 2}~ ΐ Ι ^{δ 4 +} 9 δ^{δ 6 +} ~ ] ^{y { x )}> ^{( 2}·³·^{1 7)}

^ = τ Μ ^ - £ ^ Η * * > · ^{( 2}·³·^{1 8)} For any differential operator, w e may use the expansions (2.3.13),

(2.3.16), or (2.3.18) to develop a difference relation. T o determine the

order of the truncation error, we must express the various differences in terms of the function and its derivatives, again using the Taylor series. F r o m (2.3.9) we find

(2.3.19)

(2.3.20)

(2.3.21)

(2.3.22)

(2.3.23) Powers of the difference operator are then

Similarly,

and

Using a high order approximation reduces the truncation error associated with the substitution. A s an example, consider approximating the first derivative of a function with the expression

(2.3.24)

(2.3.25)

(2.3.26) Consequently the approximation has a truncation error of order h³ whereas the simple expression

F r o m (2.3.21) we have

has a truncation error of order h.

2.4 A N A L Y T I C S O L U T I O N O F D I F F E R E N C E E Q U A T I O N S 63 It should be noted that the simple form (2.3.26) is merely a two-point formula, viz.,

T h e more complicated expression (2.3.28) requires more effort in computation and furthermore is a third order difference relation. In the next section we shall discuss the properties of solutions of difference equations and indicate additional difficulties encountered with using high order approximations.

T h u s far we have only considered approximating rather elementary differential expressions. I n later chapters we shall introduce a method of constructing appropriate difference equations for more involved equations which arise in the multigroup diffusion theory and transport theory.

2.4 Analytic Solution of Difference Equations

In order to analyze the truncation error in replacing differential equations with difference equations, w e now consider the analytic solution of difference equations. Generally speaking, the methods of solution parallel the techniques used in the differential calculus. Furthermore, the solutions of an inhomogeneous difference equation can be found by specializing the corresponding solutions of the homogeneous equation.

Except for a few simple forms of a difference equation, it is impossible to obtain a solution in closed form. Recourse is then made to a series evaluation with an attendent question of convergence. Once a problem has achieved such a degree of complexity, it is usually more satisfactory to consider numerical solutions. I n the next chapter we shall take up methods of numerically solving finite difference equations.

T h e order of a difference equation is the number of intervals separating the largest and smallest arguments of the dependent variable. T h e order is not necessarily the highest power of a difference operator. For example,

y(x + h) — y(x) dy

h ^ S * ' ^(2.3.27) whereas (2.3.24) is a four-point formula

2y(x + 3A) - 9y(x + 2h) + 18;y(x + h) - l\y(x)

6h (2.3.28)

A²y(x) — Vy(x) = y(x + 2h) - 2y(x + h) + y(x - h) = 0 (2.4.1)

and

y{x + 3Λ) - 2xy(x + h) + V(JC - h) = 3* (2.4.2) are equations of third and fourth order respectively.

A linear difference equation is one in which no products of the depend­

ent variable with itself or any of its differences appear. T h u s E q . (2.4.1) and (2.4.2) are linear, whereas

x²y(x) A*y(x) + V²v²( x ) = 0 (2.4.3) is nonlinear.

A difference equation is homogeneous if all non-zero terms involve the dependent variable; otherwise it is inhomogeneous. Equation (2.4.1) and (2.4.3) are homogeneous, but Eq. (2.4.2) is inhomogeneous. Further­

more, a difference equation may have coefficients which are constant or functions of the independent variable. I f the coefficients are constant, we speak of a difference equation with constant coefficients, otherwise we speak of an equation with variable coefficients.

2.4.1 T H E F I R S T O R D E R D I F F E R E N C E E Q U A T I O N

A first order difference equation may be simply written as

y ^ - y , = hft (2.4.4)

for which the solution is

Λ = Λ +

* Χ Λ ·

^(2-4.5)

T h e solution is given by the homogeneous solution y0 plus the particular solution, due to the inhomogeneous term/;· . A slightly more complicated first order equation would be

4vi = y*+i - yj = ^ajy> · (2.4.6)

T h e solution is obviously

Λ = Λ Π 0 . + ^a* ) . (2-4.7)

J-l

Π

J-l

• k=0

= > O e x p[ 2)^{l n}( l + « * ) ] • (2.4.8)

T h e second form of the solution corresponds to the analytic solution of the differential analog to (2.4.6).

2.4 A N A L Y T I C S O L U T I O N O F D I F F E R E N C E E Q U A T I O N S 65 As might be inferred from the above examples, all difference equations can be treated as recurrence relations. T h e approach is generally not too useful for more difficult equations, however; consequently we seek other more powerful techniques.

2 . 4 . 2 H O M O G E N E O U S D I F F E R E N C E E Q U A T I O N S W I T H C O N S T A N T C O E F F I C I E N T S

Homogeneous equations with constant coefficients can be treated by methods precisely the same as those used in differential equations. T h e general form of the difference equations would be

X = 0 . ( 2 . 4 . 9 ) 1=0

W e use the trial solution

yk = β* = e^mk. ( 2 . 4 . 1 0 ) A

Inserting ( 2 . 4 . 1 0 ) in ( 2 . 4 . 9 ) , we find nontrivial solutions exist if β is a root of the polynomial

Σ(α,)β*-<=0.

( 2 . 4 . 1 1 )

Thus, if βί are the roots of ( 2 . 4 . 1 1 )

yk=X = X b, exp [* In ft]. ( 2 . 4 . 1 2 )

i=i ; = i

T h e η coefficients bj are determined by suitable boundary conditions.

Several examples will illustrate the procedure. T h e solution of

or

( 2 . 4 . 1 3 )

( 2 . 4 . 1 4 )

is found by letting yk be given the first form of ( 2 . 4 . 1 0 ) . T h u s we have

( 2 . 4 . 1 5 )

( 2 . 4 . 1 6 )

( 2 . 4 . 1 7 )

i.e.,

T h e solution is then

β = 2 ± V2.

yk = bx{2 + V2f + 62(2 - Vlf.

y M - fyk+i + 2yk = 0, A*yk- 2Ayk-yk=0,

3β^ηι

at β = β0 if m < j - 1. Hence, — w)!] j S ^ , i.e., — m)!] β%, must be a solution of (2.4.21). I n general, therefore,

ν*=βϊ%ψ

(2-4-25)

i=0

must be a solution of (2.4.11). T h e totality of solutions consists of the roots β1, β2, ... , βη_8 plus the s independent solutions in (2.4.25).

A n example illustrates the matter. W e desire to solve

yk+2 - 4yk+1 + 4yk = 0. (2.4.26) As a second example, we note that

yk+2 - fyM + Uyk- 6y*-i = 0 (2.4.18)

has nontrivial solutions determined by the equation

β³- 6 β²+ 1 1 ) 8 - 6 = 0 . (2.4.19) Thus, the solution is

yk = cx + c22^k + c33^k. (2.4.20) As with differential equations, if one or more roots of the polynomial

Eq. (2.4.11) are repeated, then a special form of the solution results.

Suppose a root, say β0 , is repeated s times, the other roots being simple.

L e t the difference operator be represented by Θ, i.e., let

(2.4.21)

(2.4.22) where k is the independent variable. T h e n ,

where βί are the various roots of Eq. (2.4.21). Further, we note that

(2.4.23)

(2.4.24) But according to Eq. (2.4.22)

2.4 A N A L Y T I C S O L U T I O N O F D I F F E R E N C E E Q U A T I O N S 67 T h e solutions are

T h e complete solution is then

(2.4.27) Several special forms of solutions for second order equations are worth noting. I f the equation is

(2.4.28)

(2.4.29)

(2.4.30) (2.4.31)

(2.4.32)

(2.4.33)

(2.4.34) then the substitution of yk

I f b > 2, a more convenient way of writing Eq. (2.4.29) is

or

T h e solution is then

which can also be written

Another special form may be derived for the equation

for b² < 4a. I n this case, the roots are complex conjugate, and a solution is of the form

(2.4.35) where t = V— 1, ρ is the magnitude of the roots βί of Eq. (2.4.21), and θ is the arctan of the ratio of imaginary part to real part of the roots.

Alternatively

yk = p^k[c2 sin dk + c4 cos ΘΚ\. (2.4.36)

2.4.3 I N H O M O G E N E O U S D I F F E R E N C E E Q U A T I O N S

T h e solution y(j) of an inhomogeneous difference equation consists of a solution to the homogeneous equation, say yH(j)y P^{m s a} particular

solution, say yP(j), of the inhomogeneous part. I n certain cases, the particular solution may be found by the method of undetermined coefficients. W e shall not discuss this method because it is identical to the method of undetermined coefficients used with differential equations.

A more general procedure is the method of variation of parameters, which applies to difference equations with constant or variable coefficients.

L e t the difference equation be of the form

%U) = Σ P*<j)y(j + * ) = S(j), (2.4.37) fc=0

where pn(j) = 1, p0(j) Φ 0, and all other pk(j) are arbitrary. W e assume we have η linearly independent solutions for the homogeneous form of Eq. (2.4.37). Denote the homogeneous solution

= Χ^ , Ο ' ) · ( 2 A 3 8 )

i = l W e assume a particular solution of the form

3V(i)=X«,0')3\0'). (2-4.39)

1=1

In order to find the a{(j) w e need η simultaneous equations relating the

ai{])' T o accomplish this end, w e consider E q . (2.4.39) at the point j + 1; w e have

yPU + 1) = X <j + \)ylj + 1) = χ alj)yi{j + 1) + Χ ^ < ϋ> < 0 ' + i ) . i=l i = l i = l

(2.4.40) since

1) =αΜ)+ΔαΜ). (2.4.41) In order that the form of the solution be similar to E q . (2.4.39), w e try

the condition

X ^ 0 M ; + 1 ) = 0 (2.4.42)

1=1

to see where it leads. Continuing in this manner, w e generate a series of solutions of the form

yP{j + m) = X alj)yi{j + m\ (2.4.43) i = l

2.4 A N A L Y T I C S O L U T I O N O F D I F F E R E N C E E Q U A T I O N S 69

and assumed conditions

J Aali)yi{j + m) = 0, 1 < m < η - 1. (2.4.44) i=l

A t the last point w e then find

yPU + n)=j? aij)ylj + n) + £ Aatffrjj + n). (2.4.45)

W e now insert Eqs. (2.4.39), (2.4.40), etc., and Eq. (2.4.45) in the original difference equation (2.4.37). W e use the assumption that terms such as (2.4.44), vanish, and further, w e use the fact that they{( j ) satisfy the homogeneous equation. T h e result of the algebra leads to the condition

%AaMyti + n) = S{j). (2.4.46)

T h e set of assumptions [Eqs. (2.4.44) and (2.4.46)] represent a set of simultaneous equations

yiU + + + 1)^2 + . . . +yn(j + W*n = o, yx(j + 2)Aax + y2(j + 2)Aa2 + ... + yn(j + 2)Aan = 0,

(2.4.47) yx{j + n)Aax + y2(j + n)Aa2 + ... + yn(j + n)Aan = S(j).

T h e set of Eqs. (2.4.47) always has a unique solution since the deter­

minant of the coefficients, yt{j + m), cannot vanish since the yi are linearly independent by hypothesis. Once the Aai are found, the ai

are found by summation. N o t e that if S(j) = 0, then the only solution of Eq. (2.4.47) is the trivial solution as indeed it must be.

T h e procedure is illustrated by the equation

SM./) = yU + i ) - 2y(j) + y(j - 1) = /· (2.4.48) T h e homogeneous solution is

ynij) = cx + c2j. (2.4.49)

Hence w e have

yiU)

= ι

yti)=j. (2.4.50)

from which we find

= _ o' - d

«i(/) = «i(0) - 21'i=6 ² = - ^(J~Y¹- j{2) - 1), (2.4.54)

where y0 is a constant of integration. T h e particular solution is taken to be

ypU) = <i) f t [1 + 4 0 b O · (2.4.58) T h e usual condition yields

y o l l t i +c(i)]MJ) = KJ) T h e particular solution is then

yp(j) = <*i0) + (2.4.51) which leads to the equation

Aai(j) + jAa2(j) = 0,

^ i 0 ' ) + 0' + l) ^ 0' ) = i - (2.4.52) W e find that

Λ«ιΟ') = - 7². ^ 2 θ) = 7 (2.4.53)

the constant βχ(0) being taken as zero, since it may be included in the homogeneous solution anyway. Likewise,

(2.4.55) T h e complete solution is then

A s the second example we consider the most general first order equation

(2.4.56) T h e homogeneous solution by Eq. (2.4.7) is

(2.4.57)

2.5 P A R T I A L D I F F E R E N C E E Q U A T I O N S 71

or

* 0 ) = V " S W Π t¹ + ^ O l -¹. (2-4.59)

JO /c=0 t=0

the constant of " integration'' being incorporated in the homogeneous solution. T h e complete solution is thus

y(J) = Λ Π [1 + + ft Π + Φ»)] X m Π [1 + Φ")]"¹· (2-4.60) n=0 n=0 fc=0 i=0

Although the method of variation of parameters is quite general, its ability is limited to those problems for which the homogeneous solution can be found. Unfortunately the number of problems for which analytic forms of the homogeneous solution can be found is limited. Indeed, there is no general technique for solving difference equations with variable coefficients, save for rather special cases. For intractable problems approximate methods are necessary. T h e bulk of the remainder of the text is specifically devoted to methods of finding approximate solutions to difference equations.

2.5 Partial Difference Equations

Partial difference equations arise from problems involving two or more independent variables. Such equations are solved by methods similar to those used in solving partial differential equations. T h e method of separation of variables will be used in this section. W e shall consider two simple problems: an initial value problem and a boundary value problem.

For the first example consider the heat flow equation

dT(x,t) _ a²T ( x, o ) { 2 5 A

dt dx*

where Τ is the temperature. One possible difference approximation to Eq. (2.5.1) is

A T » ^T (2.5.2)

with ht as the time increment and hx the spatial increment, both assumed constant. T h e solution is desired in the semi-infinite strip 0 < χ < a.

0 < t < o o , as in Fig. 2.5.1. T h e difference equation is applied at the mesh points (k, /), that is, we desire T(k, I) = Tktl. W e shall assume the boundary conditions in the form

To.i = TK%1 = 0, all /, (2.5.3a)

and initial conditions of the form

Γ*.ο=Λ,

0 < A < # . (2.5.3b) /

1 1 1 1 t t t t

χ = 0 χ = a

F I G . 2.5.1. Mesh for the solution of the time dependent heat-flow problem.

A s in the treatment of differential equations, we assume the function Τ separable in the form

Tktl = RkSl, (2.5.4)

where Rk is a function of k only, and St a function of / only. Substituting the trial solution in the original equation, we find

htS, h*Rk ' ^{ ^

t l χ κ

where a² is an advisedly chosen separation constant and independent of k and /. T h e resulting ordinary equations in S and R may be solved by the methods of the previous section. W e find

Rk - an sin η = 1, 2,... Κ - 1, (2.5.6)

= S0e x p [ / l n ( l - <*%)], (2.5.7)

2.5 P A R T I A L D I F F E R E N C E E Q U A T I O N S 73 where S0 is an integration constant and

^ - έ - Ι ¹ - " · - ^ · ^{( 2}·⁵·^{8 )}

χ

T h e complete solution is

Tkd = Τ an sin ^ exp [/ ln (1 - <*%)], (2.5.9) where the constant S0 has been absorbed in the an. T h e coefficients an

are computed from the initial condition at / = 0. T o this end, we use the relation (see problem 11),

X sin ^ sin = - y - (1 - Sn0)8mn . (2.5.10) Multiplying (2.5.9) by sin (mnk/K) and summing over k from 0 to Κ — 1,

w e have, at / = 0,

"m = - ^ ^ f k sin , (2.5.11) which completes the solution.

In connection with the evaluation of coefficients, the following relations are useful:

cos — c o s — - = (Λ „ ^ (²·⁵·^{1 2}) m = η,

0, m Φ η, m — η even,

1, m Φ η, m — η odd,

Ό, m = η,

ο, τη Φ η> τη — η even, π(η — τη)

C Ot 2Κ +

7r(w + w ) ¹

cot ^ ,m φη,τη η

^ s , n — c o s — = j w (n + m) (2.5.13)

Κ

cos — — cos —γ- = — 8WW1(1 + δη 0) , (2.5.14)

. 2nkm . 2nkn Κ /0

2 , sin sin — — = — 6nm(\ - δη 0) , (2.5.15)

/c=0

2/ ^{c os} — g - ^{s in =} °- (2.5.16)

Ar=0

A s an example of a boundary value problem, we consider the Laplace difference equation

8Vw ι ^hi _ ο hi ^ ¹⁰ h²

y

(2.5.17) in the interval 0 < χ < α, O ^ v ^ i . W e take a uniform spacing Λχ = a/K and Ay = ft//. T h e boundary conditions are

Φθ3 = Φκό = 0,

= ο,

0fcO = gk -

^(2.5.18)

T h e mesh for this problem is shown in Fig. 2.5.2. T h e dependent variable may be separated into two functions, say R k and Sj t W e assume then,

Φkj = RkSj· (2.5.19)

O H

- J

\ t

= 0 k -= 0 k =

0 a

F I G. 2.5.2. Mesh for the solution of the Laplace difference equation.

Substituting into (2.5.17), we have

^fc+l 2 R k + Rje-l Sj- ^{2S, +}

Rk^hl

κ χ

(2.5.20) F r o m the boundary conditions, w e solve the above equations to find

flfc = s i n - ^ , » = 1,2, 1, S, = a „ s i n h ^ „ ( / - y ) ,

(2.5.21) (2.5.22)

2.6 C O N V E R G E N C E O F D I F F E R E N C E S O L U T I O N S 75

where

« ; = i^ r( l - « > s ^ ) , (2.5.23)

X

and

2/2

βη = cosh-* ( l + ^ψ-) . (2.5.24)

T h e complete solution is

0 w = Χ ^βη sinh 0n( / - /) s i n⁷^ - , (2.5.25) where

* = κ AM 1f t s i n ,x · ^{( 2}·⁵·^{2 6)} T h e resemblance to the analytic solution to Laplace's equation is noted.

2.6 Convergence of Difference Solutions

T h e convergence of the solutions of difference equations to the solu­

tions of the corresponding differential equation will be examined in this section for several simple cases. T h e differential equation is the limit of zero mesh spacing for the difference equation. T h e subject has been studied extensively, and the present state of the theory goes well beyond our present purposes. W e shall be content here with illustrating a partic­

ular method of studying convergence. T h e results will be indicative of the hazards implicit in difference approximation to differential equations.

For a simple problem, consider the ordinary differential equation

I = -«y. (2.6.D T h e solution is

y = ytf-*^x- (2.6.2)

I f we approximate the derivative by a forward difference, we have

tt+i - yi = -<*^ό · (²·⁶·³)

T h e solution is merely

y,j =y0exp [ ; Ί η( 1 - a*)]. (2.6.4)

W e expand the logarithm in the form

1 η ( 1 -αλ ) = - ο Λ - ⁽ ^ - ⁽ - ^ . (2.6.5)

I f we neglect terms beyond A², we have

y, = y0 exp [-jcA - ^ψ-] . (2.6.6)

W e identify j'A as ^ and then find

y* = Jo exp [-otXj ( l + y ) ] · (2.6.7) T h e error is thus of the form

{ah)

[ — ι - ψ - ] . (2·6·8) exp ^OLX

i.e., proportional to A. T h e smaller we take A, the smaller the error.

T h e result is consistent with the discussion of Section 2.3 where we noted that the first forward difference approximates the derivative to order A. N o t e that as A approaches 0, the difference solution converges to the differential solution.

W e might be tempted to use a more accurate difference approximation to reduce the truncation error. For instance, the approximation

Έ = ^y n — + <**"> < 2· 6· 9)

is of higher order in the truncation error. Approximating the first derivative with the difference expression in Eq. (2.6.9), w e have

= - 2 A a y , . (2.6.10) T h e solutions are of the form

yj = ^ci(V\ + a 2/ * 2 - ochy + c2(-V\ + «²A² - ochy. (2.6.11) Using the binomial expansion for the radicals, we have

V\ + <*²A² = 1 + i<*²A² - i<*⁴A⁴ + .... (2.6.12)

2.6 C O N V E R G E N C E O F D I F F E R E N C E S O L U T I O N S 77

I f terms beyond A⁴ are neglected, the first term becomes

yiU) = (ι - ο* + - y - - i « η , ( 2 . 6 . 1 3 )

= βχ ρ/ [ - α Λ + - ^ ] , ( 2 . 6 . 1 4 )

where w e have truncated the logarithmic expansion to obtain ( 2 . 6 . 1 4 ) . N o t e that the function yx indeed approximates the analytic solution to order h² and is more accurate than Eq. ( 2 . 6. 7 ) for the same h. However, reference to the second term of Eq. ( 2 . 6 . 1 1 ) shows that the second independent solution is of the form

yjj) = ( - y( V l + c W + c*)'. ( 2 . 6 . 1 5 ) For all h > 0 this function is ever increasing in magnitude and oscillatory

in sign. Hence, this second form of the solution would completely domi­

nate the decaying exponential of the first term, and an entirely erroneous solution would result. It might be argued that the coefficient c2 = 0, and hence this difficulty does not arise. However, for numerical computa­

tions any round-off might introduce the second independent solution and continued calculation would result in a complete degeneration of the problem.

T h e basic difficulty with the higher order approximation is the introduction of additional independent solutions. I n many cases various procedures may be adopted to control the behavior of such extraneous solutions, and we shall consider such problems in later chapters. For now w e merely note one of the difficulties encountered in attempting to find high order accuracy in our solution.

A s a more detailed study of convergence, we consider the boundary value problem of Section 2 . 5 . T h e solution ( 2 . 5 . 2 5 ) contained the term sin (nnk/K) . In the limit w e have

.. . rrnk . nnkhx . ηπχ

hm sin —rr = lim sin = sin , ( 2 . 6 . 1 6 )

hx-*o Κ hx->o a a ^{v 7}

which is the analytic result for this factor. For the separation constants, we have

ga-s-ea [m "-H (v)'«]- <"·">

( " . . 8 )

which is the analytic result. Further,

βη

so that

Therefore,

l i m c o s h ^ = l i m [ l + - ^ ] , (2.6.19)

TlTT

Urn βη = lim α J L = lim hy . (2.6.20)

TL7T

lim sinh /?„(/ - j) = sinh — (b - y), (2.6.21)

ΛΒ-»0 U

which is the correct result. I n other words, the difference solution does indeed converge to the differential solution. T h e truncation error in the solution appears in the separation constants and is of the form

an (difference ^{= a}n [differential [ l ~~ {~α~ί ' (2.6.22)

βη difference ⁼ βη I differential ~~ J2 i~a~) ^ ] * (2.6.23) T h e truncation error is proportional to the square of the mesh spacing since central differences were used.

For the final study w e consider the exponential term in the heat-flow example of the previous section. F r o m E q . (2.5.9) w e have

lim /In (1 - *X) = lira / In [ l - ^ ( l - cos £ ) ] . (2.6.24)

U p o n expanding the logarithm in the left-hand side of (2.6.24), w e have l i m / ( - « « * . ( 2 . 6 . 2 5 )

I f w e take t = lhty then the first term becomes

e-«lt = (2.6.26) which is the analytic solution. F r o m this w e might conclude the solution

converges. However, examination of the right-hand side of (2.6.24) shows that for large η (η ^ Κ — 1 by E q . (2.5.6)), the term

1 - c o s - ^ ^ + 2 . (2.6.27)

2.7 M A T R I X F O R M O F D I F F E R E N C E E Q U A T I O N S 79

T h e right-hand side becomes

lim/1„ [ l _ _ * * « . ] . (2.6.28) I f hx and ht are related such that ht/hx2 > ^, then the argument in

(2.6.28) is of magnitude greater than unity. T h e logarithm is positive (plus a phase factor) and hence the solution will diverge, thereby entirely failing to represent the desired solution. I n order to guarantee the con­

vergence of the approximation, w e must require

\ < \ (2-6.29) even as ht,hx-+ 0.

Requirements such as (2.6.29) are very frequent in numerical studies, and w e shall encounter such relations often in our later analysis.

2.7 Matrix Form of Difference Equations

In our later work w e shall frequently use a compact matrix notation for the simultaneous equations resulting from difference approximation to differential equations. T o illustrate the construction of the matrix form of the equation, w e n o w consider the t w o examples of partial difference equations given in Section 2.5.

W e approximate the heat-flow equation by the simple difference relation

Tk,i+i — Tktl — -j— [Tk+ltl — 2Tktl + Tk_ia], (2.7.1)

X

W e assume the boundary conditions as before, Eqs. (2.5.3a,b). Equation (2.7.1) can be factored in the form

h

Tk.i+i = -JJT [Tk+ltt — ocTktl + Tk_ltl] (2.7.2) with a = 2 — hx2/ht.

T h e various equations for a given / are then - p - ( — α Γ1 >ζ + T2fl) = Tl t l + 1,

X

h

~ ^ 2 ~ ( ^ l . Z ^{— a}^ 2 . i + ^ 3 . 0 ⁼ ^2.1+1 >

-W ( Γ * - . . , - «ZV-i.i) = r H i . (2.7.3c) W

(2.7.3a)

(2.7.3b)

T h e above set of simultaneous equations can evidently be written

- α 1 0 ... 0 0"

1 - α 1 ... 0 0

0 0 0 ... -OL 1 0 0 0 ... 1 -oc

L^V-I.ZJ

L^tf-l.l+lJ

ι,ι+ι

^2,1+1

(2.7.4) I f w e denote the matrix in (2.7.4) as A , and define the vector ψ , as

T2l

(2.7.5)

then the entire set of equations can be written

hi Α ψ( = ψί + 1. (2.7.6)

I f w e denote the starting vector as ψ0 , then E q . (2.7.6) can also be written

[(Α,/ΑΪΪΑΐ'ψο = ψ , .

N o t e that the matrix A is of tridiagonal form.

For the Laplace difference equation, w e have the relation

h² ^ h²

(2.7.7)

(2.7.8)

By simple factoring, w e have

r2tyk+lJ +Φΐ€-υ) — βΦν + ΦΜ+1 + ΦkJ-l = 0,

where r² - h²\h,² and β = 2(1 + r²).

(2.7.9)

2.7 M A T R I X F O R M O F D I F F E R E N C E E Q U A T I O N S

W e define the vector ψ^· as

Φΐ)

<l>2j

.Φκ-υ J

T h e set of equations (2.7.9) is then

Α ψ , + Ι ψί +1 + Ι ψ ^ = 0, where A is the matrix

A =

-β r² 0 ... 0 0"

r² -β r² ... 0 0

0 0 0 ... r² -β_

W e now define the extended vector ψ

Ψ =

"Ψι Ψ2

Ψ/-ι

ΓΦη

Φ21

Φΐ2 Φ22

L<^Ar-i,y-i-J

T h e set of equations (2.7.11) can then be written Β ψ = 0 ,

where

Β =

Γ Α I 0 ... 0 I A I ... 0

I A

81

(2.7.10)

(2.7.11)

(2.7.12)

(2.7.13)

(2.7.14)

Each of the elements of Β is a (Κ — 1) by (Κ — 1) square submatrix.

A matrix of the form of Β is a block tridiagonal matrix. T h e procedure can obviously extended to difference equations with more unknowns and also to equations with variable coefficients and/or variable mesh spacing.

References

There are several books devoted to the calculus of finite differences. Parti­

cularly exhaustive treatments are included in references 7, 2, and 3. Somewhat shorter treatments are found in 4 and 5. Reference 5 also includes many other aspects of the difference calculus. T h e application of finite differences is not limited to approximating differential equations. Other uses for the calculus are described very cogently in 4 and 6.

L Jordan, C , "Calculus of Finite Differences." Chelsea, N e w York, 1947.

2. Milne-Thompson, L. M . , " T h e Calculus of Finite Differences.'' Macmillan, London, 1933.

3. Fort, T., "Finite Differences and Difference Equations in the Real Domain.'' Oxford Univ. Press, London and N e w York, 1948.

4. Hildebrand, F. B., "Introduction to Numerical Analysis." McGraw-Hill, N e w York, 1956.

5. Hildebrand, F. B., "Methods of Applied Mathematics." Prentice-Hall, Englewood Cliffs, N e w Jersey, 1952.

6. Hamming, R. W . , "Numerical Methods for Scientists and Engineers." M c G r a w - Hill, N e w York, 1962.

Problems

1. Prove the following:

. . . _ . ah ι , ah\

(a) Δ sm{ax + b) = 2 sin — cos + b + — j , (b) Δ cos(tf# - f b) = —2 sin — sin \ax + b + — j . 2. Show that

/(* + nh) = /(*) + ηΔ/(χ) + ^{n (n} ~ ^1} j y ( * ) + ... + Af(x).

3. Show that

^ sin(x — h/2) (a) £ cos 2

(b) ^ sinh χ =

2

sin

h/2

cosh(* — h/2) 2 sinh h/2

4. Prove that the differentiation operator d/dx may be formally written

— = τ sinh ¹ -

dx h

2

and then develop expansion (2.3.17).

P R O B L E M S 8 3

5. Find the solution to the following:

(a) A*gk + Agk - gk = 0, (b) g(k + 3) - 8*(Λ) = 0.

6. Find the solutions to the eigenvalue problem + A«y* = 0, 0

g(Q) = l,

with boundary conditions

4Vo = 0; Vy* = 0.

Show the convergence to the analytic solution.

A n infinite medium, void of neutrons, has a constant scattering cross-section and a constant fission cross-section. There is no capture. T h e fission process produces 1 prompt neutron and another neutron delayed exactly τ sees.

T h e total diffusion time of any neutron is exactly τ sees. One neutron is introduced at time t = 0 into the assembly.

[a] Derive an expression for the neutron population as a function of time.

[b] T h e assembly is to be scrammed when more than 10¹⁵ neutrons are present. W h e n should scram occur ?

8. Approximate the wave equation 3*φ dt² with the difference approximation

θ²φ

8²φ

Solve the difference equation and compare the solutions with the analytic solution to the differential equation. Under what circumstances will the difference equation converge ?

9. The determinant, ΔΝ, of the Ν by Ν matrix - « 1 0 ·

1 - a 1 0

arises from applying the operator δ². By assuming A0 = 1, show the recurrence relation

is valid.

AN = - AN_2

10. Derive a difference approximation to the Laplacian operator in r, ζ coordi­

nates which is accurate to order h\ and h\.

11. By noting that the series Σ^ο* e^inkn/K is geometric, prove

(a)

(b)

[The other orthogonality relations given in Section 2.5 may be developed in a similar manner].

12. The purpose of this problem is to illustrate the construction of a difference equation by physical rather than mathematical means. Heat is transferred from a cylindrical fuel element of radius R() encased in cladding of uniform thickness a, centered in a cylindrical coolant channel of radius RK . The coolant channel is assumed insulated at the outer boundary; the coolant mass flow rate and inlet temperature are assumed fixed. Calculate the tem­

perature of the fuel, cladding, and coolant as functions of r and z> the radial and axial coordinates. Neglect any heating in the coolant or cladding as a result of neutron or photon irradiation. Assume thermal conductivities are independent of temperature in all materials; neglect axial heat flow; neglect any thermal resistance between cladding and fuel. Assume the heat source in the fuel depends only on axial position and that the coolant is completely turbulent.

(a) Write down the usual energy conservation equations and boundary conditions needed to solve the problem. Let η be the heat transfer coefficient between the cladding and the coolant, and k^f the fuel thermal conductivity.

(b) Solve the resulting equations analytically for the temperature as a function of r and ζ in the fuel, cladding, and coolant.

(c) Divide the assembly into a number of annuli such that the radius tj of the 7th one is given by r, = jh. By performing energy balances directly on zones j — 1, j\ and j + 1, show that

-^i-

(Γ<_ι - T ) + ^) - ( Τί +1 - T,) =

S (zW) - rU,

where T} is temperature at rt and S(z) is the heat source density. What happens at j = 1, i.e., how is the temperature of the centerline to be deter­

mined ? Show that the above result is correct to first order.

(d) By using the cross-sectional area of flow at the average radius of an interval, show that

- ^ - ( T , - ! - T^r^ + r,) + - £ - ( Τ ,+1 - Τ, ) (Γ| + = S(z)(r) - rJ.J.

Eliminate r^x, rj, and ri+1 in terms of j and h. Show that the truncation error is of order h.

(e) As the next improvement, for the source use the volume between tj — h/2 and r} + h/2f and show that the heat flow equation is

- £ - [ ( T f- i - T^r^ + r,) + ( T ,+l - T,)(r, + r ,+ 1) ] = S ( * ) ( r J+i - rj.j).

P R O B L E M S 85 Eliminate the various r} in terms of j and h. Show that the truncation error is of order h².

(f) Use central differences in the result of (a) and show that the result agrees with that (e).

(g) Determine the center temperature in (e) in terms of Tx and q(z) by evaluating limr-*o [(\/r)(dT/dr)]. Hint: Expand (dT/dr)\r=0 in a Taylor series about r = 0. Show that

(h) H o w should the cladding be treated ?

(i) H o w should the equations be generalized if the coefficients depend on r, zy and T?