Finite Difference Methods

Research is to see what everybody else has seen, and think what nobody has thought. Albert SzentGyorgyi

I. Introduction Analytical methods may fail if: 1. The PDE is not linear and cant be linearized without seriously affecting the result. 2. The solution region is complex. 3. The boundary conditions are of mixed types. 4. The boundary conditions are time-dependent. 5. The medium is inhomogeneous or anisotropic. The finite difference method (FDM) was first developed by A. Thom* in the 1920s under the title the method of square to solve nonlinear hydrodynamic equations. * A. Thom an C. J. Apelt, Field Computations in Engineering and Physics. London: D. Van Nostrand, 1961. The finite difference techniques are based upon the approximations that permit replacing differential equations by finite difference equations. These finite difference approximations are algebraic in form, and the solutions are related to grid points. Thus, a finite difference solution basically involves three steps: 1. Dividing the solution into grids of nodes. 2. Approximating the given Fig. 1 Common two-dimensional grid differential equation by finite patterns difference equivalence that relates the solutions to grid points. 3. Solving the difference equations subject to the prescribed boundary conditions and/or initial conditions. II. Finite Difference Scheme Differential equations equations estimating derivatives numerically finite difference

Given a function f(x) shown in Fig. 2, we can approximate its derivative, slope or tangent at P by the slope of the arcs PB, PA, or AB, for obtaining the forwarddifference, backward-difference, and central-difference formulas respectively.

Finite-Difference-Method-for-PDE-1

forward-difference formula f ( xo + x ) f ( xo ) f ' ( xo ) x

backward-difference formula f ( x o ) f ( x o x ) f ' ( xo ) x central-difference formula f ( xo + x ) f ( xo x ) 2x

f ' ( xo )

Fig. 2 Estimates for the derivative of f(x) at P by using forward, backward, and central differences.

The approach used for obtaining above finite difference equations is Taylors series: 1 1 f ( xo + x ) = f ( xo ) + xf ' ( xo ) + ( x ) 2 f '' ( xo ) + ( x ) 3 f ''' ( xo ) + O ( x ) 4 , (1) 2! 3! and 1 1 f ( xo x ) = f ( xo ) xf ' ( xo ) + ( x ) 2 f '' ( xo ) ( x ) 3 f ''' ( xo ) + O ( x ) 4 , (2) 2! 3! where O(x)4 is the error introduced by truncating the series. To subtract (1) by (2), we can obtain f ( xo + x ) f ( xo x ) = 2xf ' ( xo ) + O ( x ) 3 , which could be re-written as f ( x o + x ) f ( x o x ) f ' ( xo ) + O ( x ) 2 , i.e. the central-difference formula. Note 2x 2 that the O(x) means the truncation error is the order of (x)2 for the centraldifference. The forward-difference and backward-difference formulas could be obtained by re-arranging (1) and (2) respectively, and we have f ( xo + x ) f ( xo ) f ' ( xo ) + O ( x ) , for forward difference, x and f ( x o ) f ( x o x ) f ' ( xo ) + O ( x ) , for backward difference. We can find the x truncation errors of these two formulas are of order x. Upon adding (1) and (2), f ( x o + x ) f ( x o x ) = 2 f ( x o ) + ( x ) 2 f '' ( x o ) + O ( x ) 4 , and we have f ( x o + x ) 2 f ( x o ) + f ( x o x ) f '' ( x o ) + O ( x ) 2 . ( x ) 2 Higher order finite difference approximations can be obtained by taking more terms in Taylor series expansion.

Finite-Difference-Method-for-PDE-2

To apply the difference method to find the solution of a function (x,t), we divide the solution region in x-t plane into equal retangles or meshes of sides x and t. The derivatives of at the (i,j)th node are shown in the table, where x = ix t = jt .

Finite-Difference-Method-for-PDE-3

Finite Differencing of Parabolic PDEs

Consider a simple example of a parabolic (or diffusion) partial differential equation with one spatial independent variable 2 k (3) = 2 , t x where k is a constant. The equivalent finite deference approximation is (i, j + 1) (i, j ) (i 1, j ) 2 (i, j ) + (i + 1, j ) . (4) k = t ( x ) 2 where x=ix, i=1,2,3,,n, t=jt, j=1,2,. In (4), we use the forward difference formula for the derivative with respective to t and central difference formula for the with respect to x. If we let t r= , (5) k ( x ) 2 Eq. (4) could be written as (i, j + 1) = r (i + 1, j ) + (1 2r ) (i, j ) + r (i 1, j ) . (6)

This explicit formula can be used to compute (x,t+t) explicitly in terms of (x,t). Thus the values of along the first time row (see Fig.3), t=t, can be calculated in terms of the boundary and initial conditions, then the values of along the second row, t=2t, are calculated in terms of the first time row, and so on.
Fig. 3 Finite difference mesh for two independent variable x and t.

A graphic way of describing the difference equation (6) is through the computational molecule of fig. 4, where the square is used to represent the grid point where is presumed known and a circle where is unknown.

Fig. 4 Computational molecule for parabolic PDE: (a) for 0 < r < 1/2 (b) r = 1/2.

Finite-Difference-Method-for-PDE-4

Stability Analysis of the explicit algorithm

Reducing the mesh size could increase accuracy, but the mesh size could not be infinitesimal. Decreasing the truncation error by using a finer mesh may result in increasing the round-off error due to the increased number of arithmetic operations. A point is reached where minimum total error occurs for any particular algorithm using any given word length. This is illustrated in Fig. 5.

Fig. 5 Error as a function of the mesh size.

The concern about accuracy leads us to a question whether the finite difference solution can grow unbounded, or a property termed the instability of the difference scheme. A numerical algorithm is said to be stable if a small error at any stage produces a smaller cumulative error. Otherwise, it is unstable. To determine a whether a finite difference scheme is stable, we define an error, n, which occurs at time step n, assuming there is one independent variable. We define the amplification of this error at time step n+1 as (7) n+1=gn, where g is the amplification factor. In more complicated situation, we have more independent variables, and (7) becomes (8) []n+1=[G][]n, where [G] is amplification matrix. For the stability of the finite difference scheme, it is required that Eq. (7) satisfy (9) |n+1| |n|, or (10) |g| 1. Forthecasein(8),thedeterminantof[G]mustvanish,i.e. Det[G]=0. (11)

One useful method and simple method of finding a stability criterion for a finite difference scheme is to construct a Fourier analysis of the difference equation and thereby derive the amplification factor. The technique is known as von Neumanns method. The stability condition is called the von Neumann condition. (ref. To J. W. Thomas, Numerical Partial
Differential Equations Finite Difference Methods, Springer-Verlag New York, 1995)

Finite-Difference-Method-for-PDE-5

Considering the explicit scheme of Eq. (6): in +1 = (1 2r ) in + r ( in+1 + in1 ) ,

(12)

where r=t/k(x)2. We have changed our notation so that we can use j = 1 in the Fourier series. Suppose the solution is in = An (t )e jix , 0 x 1, (13) where is the wave number. Since the differential equation is linear, we need consider only one Fourier mode, i.e. in = An (t )e jix . (14) Substituting (14) into (12) gives An +1e jix = (1 2r ) An e jix + r ( e jx + e jx ) An e jix or An +1 = An [1 2r + 2r cosx ] .

(15) (16)

Hence the amplification factor is An +1 x g = n = 1 2r + 2r cos x = 1 4r sin 2 . (17) A 2 In order to satisfy (10) x 1 4r sin 2 1. (18) 2 Since this condition must hold for every wave number , we take the maximum value of the sine function so that 1 4r 1 and r 0 or r and r 0. Of course, r = 0 implies t = 0, which is impractical. Thus we have 0 < r 1/2. In order to ensure a stable solution or reduce errors, care must be exercised in selecting the value of r in (5) and (6). If we choose r=1/2, then (12) becomes 1 in +1 = ( in+1 + in1 ) , (19) 2 as shown in Fig. 4(b).

Finite-Difference-Method-for-PDE-6

Alternative Schemes for parabolic PDEs

Although the formula is simple to implement, its computation is slow. An implicit formula, proposed by Crank and Nicholson in 1974, is valid for all finite values of r. 2 in Eq. (3) by the average of the central difference formulas on the jth We replace x 2 and (j+1)th time rows such that (i, j + 1) (i, j ) 1 (i 1, j ) 2 (i, j ) + (i + 1, j ) = [ k 2 ( x ) 2 t , (20) (i 1, j + 1) 2 (i, j + 1) + (i + 1, j + 1) ] + ( x ) 2 which can be re-written as r (i 1, j + 1) + 2(1 + r ) (i, j + 1) r (i + 1, j + 1) = r (i 1, j ) + 2(1 r ) (i, j ) + r (i + 1, j ) . (21)

Where r is defined in (5). The left hand side of (21) consists of three unknown values of . This is illustrated in the computational molecule of Fig. 5(a). Thus if there are n nodes along each time row, then for j=0, applying (21) to nodes i=1,2,,n results in n Fig. 5 Computational molecule for Cranksimultaneous equations with n Nicholson method: (a) for finite r and (b) r = 1. unknown values of and known initial and boundary values of . Similarly, for j=1, we obtain n simultaneous equations for n unknown values of in terms of the known values at j=0, and so on. The combination of accuracy and unconditional stability allows the use of a much larger time step with Crank-Nicholson method than is possible with the explicit formula. Although the method is valid for all finite values of r, a convenient choice of r=1 reduces (21) to (i 1, j + 1) + 4 (i, j + 1) (i + 1, j + 1) = (i 1, j ) + (i + 1, j ) , with the computational molecule of Fig. 5(b). <HW> Use the von Neumann approach to determine the stability condition of Eq. (21). (22)

Finite-Difference-Method-for-PDE-7

Finite-Difference-Method-for-PDE-8

[Example] Solve the diffusion equation 2 0x1 = x 2 t subject to the boundary conditions (0,t) = 0, (1,t) = 0, t > 0 and initial condition (x,0) = 100. Solution This problem may be regarded as a mathematical model of the temperature distribution in a rod of length L=1m with its end in contacts with ice blocks (or held at

0) and the rod initially at 100. With the physical interpretation, out problem is finding the internal temperature as a function of position and time. We will solve this problem using both explicit and implicit methods. (a)Analytical solution ( x, t ) =

sin(nx ) exp( n 400 n

2

t) ,

n=2k+1

k =0

t x 0.00 0.10 0.20 0.30 0.40 0.50 ------------------------------------------------0.0000 50.00 100.00 100.00 100.00 100.00 100.00 0.0050 0.00 68.27 95.45 99.73 99.99 100.00 0.0100 0.00 52.05 84.27 96.61 99.53 99.92 0.0150 0.00 43.63 75.18 91.67 97.85 99.22 0.0200 0.00 38.29 68.26 86.59 95.18 97.52 0.0250 0.00 34.52 62.86 81.85 91.91 94.93 0.0300 0.00 31.67 58.47 77.51 88.32 91.75 0.0350 0.00 29.39 54.78 73.50 84.61 88.24 0.0400 0.00 27.50 51.58 69.78 80.88 84.58 0.0450 0.00 25.87 48.74 66.31 77.21 80.88 0.0500 0.00 24.42 46.16 63.04 73.63 77.23 0.0550 0.00 23.12 43.79 59.96 70.18 73.67 0.0600 0.00 21.93 41.59 57.04 66.86 70.22 0.0650 0.00 20.82 39.53 54.27 63.68 66.90 0.0700 0.00 19.79 37.59 51.65 60.63 63.72 0.0750 0.00 18.81 35.75 49.15 57.73 60.68 0.0800 0.00 17.89 34.01 46.78 54.96 57.78 0.0850 0.00 17.02 32.37 44.52 52.32 55.00 0.0900 0.00 16.20 30.80 42.38 49.81 52.36 0.0950 0.00 15.41 29.31 40.34 47.41 49.85 0.1000 0.00 14.67 27.90 38.39 45.13 47.45

(b) Explicit Method For simplicity, let us choose x=0.1, r=1/2 so that r ( x ) 2 t = = 0.005 k since k=1. We need the solution for only 0 x 0.5 due to the fact that the problem is symmetric with respect to x=0.5. Notice that the value of (0,0) and (1,0) are taken as the average of 0 and 100.

Finite-Difference-Method-for-PDE-9

t x 0.00 0.10 0.20 0.30 0.40 0.50 ------------------------------------------------0.0000 50.00 100.00 100.00 100.00 100.00 100.00 0.0050 0.00 75.00 100.00 100.00 100.00 100.00 0.0100 0.00 50.00 87.50 100.00 100.00 100.00 0.0150 0.00 43.75 75.00 93.75 100.00 100.00 0.0200 0.00 37.50 68.75 87.50 96.88 100.00 0.0250 0.00 34.38 62.50 82.81 93.75 96.88 0.0300 0.00 31.25 58.59 78.12 89.84 93.75 0.0350 0.00 29.30 54.69 74.22 85.94 89.84 0.0400 0.00 27.34 51.76 70.31 82.03 85.94 0.0450 0.00 25.88 48.83 66.89 78.12 82.03 0.0500 0.00 24.41 46.39 63.48 74.46 78.12 0.0550 0.00 23.19 43.95 60.42 70.80 74.46 0.0600 0.00 21.97 41.81 57.37 67.44 70.80 0.0650 0.00 20.90 39.67 54.63 64.09 67.44 0.0700 0.00 19.84 37.77 51.88 61.04 64.09 0.0750 0.00 18.88 35.86 49.40 57.98 61.04 0.0800 0.00 17.93 34.14 46.92 55.22 57.98 0.0850 0.00 17.07 32.42 44.68 52.45 55.22 0.0900 0.00 16.21 30.88 42.44 49.95 52.45 0.0950 0.00 15.44 29.33 40.41 47.45 49.95 0.1000 0.00 14.66 27.92 38.39 45.18 47.45

(c) Implicit Method Let us choose x=0.1, r=1 so that t=0.01.

t x 0.00 0.10 0.20 0.30 0.40 0.50 ------------------------------------------------0.0000 50.00 100.00 100.00 100.00 100.00 100.00 0.0100 0.00 59.81 89.23 97.10 99.17 99.59 0.0200 0.00 40.18 71.50 88.92 95.77 97.47 0.0300 0.00 32.96 60.34 79.30 89.59 92.68 0.0400 0.00 28.33 52.97 71.30 82.31 85.95 0.0500 0.00 25.06 47.28 64.41 75.09 78.70 0.0600 0.00 22.46 42.56 58.30 68.27 71.68 0.0700 0.00 20.25 38.46 52.82 61.98 65.13 0.0800 0.00 18.32 34.82 47.87 56.23 59.11 0.0900 0.00 16.59 31.54 43.40 51.00 53.61 0.1000 0.00 15.03 28.59 39.34 46.24 48.62
50 40

temperature

30 20 10 0 Analytical solutions Explicit method Implicit method 0 0.1 0.2 0.3 0.4 0.5

The temperature profiles at t=0.1

Finite-Difference-Method-for-PDE-10

Finite Differencing of Hyperbolic PDEs

Consider a simple example of a hyperbolic partial differential (or wave) equation with one spatial independent variable 2 2 2 (23) = 2 , u x 2 t where u is the speed of the wave. The equivalent finite deference approximation is (i 1, j ) 2 (i, j ) + (i + 1, j ) (i, j + 1) 2 (i, j ) + (i, j 1) u2 . (24) = ( x ) 2 ( t ) 2 where x=ix, i=1,2,3,,n, t=jt, j=1,2,. In (24), we use the central difference formula for the derivatives with respective to t as well as with respect to x. If we let 2 ut r= (25) , x Eq. (24) could be written as (i, j + 1) = 2(1 r ) (i, j ) + r[ (i + 1, j ) + (i 1, j )] (i, j 1) . (26) If we choose r = 1, Eq. (26) becomes (i, j + 1) = (i + 1, j ) + (i 1, j ) (i, j 1) . (27)

(28) (26) (29)

(30) (30) (31) (31)

(32)

Finite-Difference-Method-for-PDE-11

[Example]

Solution

(33) (26) (34)

(34) (35) (34) I (33) (35) I

Finite-Difference-Method-for-PDE-12

Table I

Figure I

Finite-Difference-Method-for-PDE-13

Finite-Difference-Method-for-PDE-14

References
1. "Numerical Techniques in Electromagnetics," by Matthew and Sadiku, CRC Press, Inc. (1992) 2. Applied Numerical Analysis, by C. F. Gerald and P. O. Wheatley, Addision Wesley Longman, Inc. (1997) 3. High Performance Computing, by Kelvin Dowd and C. R. Severance, OReilly and Associates, Inc. (1998)

Finite-Difference-Method-for-PDE-15