Download Partial Differential Equations - Seminar in Engineering Analysis - Lecture Notes and more Study notes Engineering Analysis in PDF only on Docsity!
Numerical Analysis of Partial Differential Equations
Solution Properties for Finite-difference Equations
A numerical solution to an ordinary or a partial-differential equation should satisfy various properties. These are listed below.
Consistency. A finite-difference expression is consistent with the differential equation if the truncation error approaches zero (ignoring roundoff) as the discrete steps (in space and time) approach zero. This is usually the case for finite-difference expressions. However, there are some algorithms, such as the DuFort-Frankel algorithm, in which the truncation error depends on the ratio (Δt/Δx)^2 , that are only conditionally consistent.
Stability. A finite-difference equation should be stable. The errors in a stable finite-difference equation will not grow without bound as the solution progresses.
Convergent. A convergent finite-difference expression tends to the exact solution of the differential equation as the grid size tends towards zero. According to the Lax Equivalence Theorem, a properly posed linear initial value problem and a finite difference approximation to it that satisfies the consistency condition will be convergent if it is stable.
Physical reality. Solutions should produce physically realistic results. Densities should be positive. Temperature changes should not violate the second law of thermodynamics. This should be true for each node in the solution. This requirement applies not only to the numerical method, but to physical models for complex flow phenomena such as turbulence, combustion, gaseous radiation, and multiphase flow.
Accuracy. There are many sources of error in numerical solutions. We have discussed truncation errors caused by the numerical approaches. Additional errors, known as iteration errors, are possible when approximate solutions to the finite-difference equations are obtained. Furthermore, inaccuracies can be introduced by poor physical models or assumptions. For example, the solution of a potential flow problem will have possible errors from truncation and iteration error. However, the assumption of potential flow can introduce errors as compared to the actual physical problem. However, these errors may be acceptable in some cases.
Finite-difference methods and stability for the diffusion equation
The notes on the entitled “Introduction to Numerical Calculus”, referred to here as “introductory notes”, we applied finite-difference and finite-element methods to a simple ordinary differential equation. The extension of the finite-difference approach used there to partial differential equations is fairly straightforward. As an example of this, consider the following differential equation, known as the diffusion equation or the heat conduction equation.
2
2
x
T
t
T
[3-1]
The quantity represents the thermal diffusivity in heat transfer and the diffusion coefficient in diffusion problems. We will call this term the diffusivity. The dependent variable T can be a general potential, although here we are using the usual symbol for temperature. Equation [3-1]
has an open boundary in time; we do not have to specify the conditions at a future time. This equation is formally classified as a parabolic partial differential equation.
We need an initial condition for the temperature at all points in the region. We can write a general initial condition as T(x,0) = T 0 (x). Similarly we can write arbitrary boundary conditions for the boundaries at xmin and xmax: T(xmin,t) = TL(t) and T(xmax,t) = TR(t). These conditions, as well as values for the geometry and the diffusivity must be specified before we can solve the differential equation.
We can construct finite-difference grids in space and time. For simplicity, we will assume constant step sizes, Δt and Δx. We define our time and space coordinates on this grid by the integers n and i, so that we can write
t n t 0 n t and xi x 0 i t [3-2]
We can obtain a solution process, known as the explicit method , if we use a forward difference for the time derivative and a central difference for the space derivatives. This method is also called the forward-time-central-space (FTCS) method to emphasize the nature of the finite- difference approximations to the time and space derivatives in the diffusion equation. The forward-time and central-space approximations are given by the following equations.
[( )]
(^12)
O x
x
T T T
x
T
O t and
t
T T
t
T
n i
n i
n i
n
i
n i
n i
n
i
[3-3]
These equations are modifications of equations [19] and [29] from the introductory notes on numerical calculus. In those notes we dealt with ordinary derivatives and needed only one subscript for the dependent grid variable. Here, in a multidimensional case, we have dependent variable, T, as a function of time and distance, T(x,t). Thus we define Tin^ = T(xi, tn). We use differences between variables like Tin, T ni+1, and Tni-1 to give differences in the x-direction that are used to compute spatial derivatives. The constant n superscript in the space derivative is an indication that we are holding time constant and varying distance. Similarly, the time derivatives are reflected in the finite-difference form by using variations in the superscript n while holding the subscript i constant.
The forward difference for the time derivative is chosen over the more accurate central difference to provide a simple algorithm for solving this problem. There are other problems with the use of the central difference expression for the time derivative that we will discuss later.
If we substitute the finite-difference expressions in [3-3] into equation [3-1], we get the following result.
[ ,( )]
2
1 1
1
O t x
x
T T T
t
Ti n Tin in in in
[3-4]
We have written the combined error term to show the order of the error for both Δx and Δt. We have not written these as separate terms since O() notation simply shows how the error depends on the power of the step size. If we neglect the truncation error, we can obtain a simple equation for the value of the temperature at the new time step.
i n T^ in Tin Tin ^ Tin
x
t
T
2 1 1
1 [3-5]
This gives Nx = 4 and Δx = 0.25. We also pick a time step, Δt = 0.01 and calculate to a maximum time of 0.2. This requires 20 time steps.
Table 3-1 below shows the results of the numerical calculations with the spatial variation in columns and the time variation in rows. There are two initial condition rows. The “t = 0” row shows the initial conditions. The “t = 0+” row shows the temperatures immediately after t = 0 when the boundary conditions are applied. This second row is the starting time point, n = 0, for the numerical calculations. These two initial condition rows illustrate the discontinuity in the T values at the two points (0,0) and (L,0) in the solution domain. This discontinuity does not affect our ability to obtain an analytical solution. We must treat it explicitly in the numerical solution.
In addition to the modified initial condition row, where we know all the values of T, we have the boundary condition columns. For all time points at x = 0 and x = L = 1, we set the table entries to the known boundary values. Although the initial and boundary conditions are particularly simple in this example, we could handle initial conditions that were arbitrary functions of x and boundary conditions that were arbitrary conditions of time by placing the appropriate results for the initial and boundary temperatures at the correct places in the table.
Once the modified initial condition row and the boundary condition columns are set, the solution for the remaining temperatures, using equation [3-7] is straightforward. We first use equation [3- 6] to compute the value of f.
( )^2
x
t
f
[3-9]
Once the value of f is known, we can apply equation [3-7] to each node, in successive rows. The calculations below are for the first time step.
( 1 2 ) 0. 16 [ 1000 0 ] 0. 68 [ 1000 ] 840
( 1 2 ) 0. 16 [ 1000 1000 ] 0. 68 [ 1000 ] 1000
( 1 2 ) 0. 16 [ 0 1000 ] 0. 68 [ 1000 ] 840
0 3
0 4
0 2
1 3
0 2
0 3
0 1
1 2
0 1
0 2
0 0
1 1
T fT T f T
T fT T f T
T fT T f T
[3-10]
These are the results shown in the “n = 1” row of the table. The numerical results are symmetric in x as would be expected from the uniform initial condition and the symmetric boundary condition. The process illustrated in equation [3-10] is used to calculate the remaining results in the table. The final two entries in the table are the exact results, calculated for the different values of x at the final time t = 0.2, and the error. The latter is the absolute value of the difference between the numerical and analytical solution.
The numerical calculation of the final values for an end time of t = 0.2 can be repeated for time steps of 0.02 and 0.04. These results are shown in the Tables 3-2 and 3-3 below. When the time step is doubled from 0.01 to 0.02 the solution results are similar. However, the error in the temperatures at the final time step has increased by about a factor of three.
When the time step is increases to 0.04, the results are completely different, and these results do not match our ideas of physical reality. We expect that the region will eventually reach a constant value of T equal to the two boundary temperatures of zero. As this diffusion process takes place, we do not expect to see any T values outside the range bounded by the initial value of 1000 and the boundary value of 0.
Numerical solution of the Diffusion Equation
- Table 3- - = 1, Δx = 0.25, Δt = 0.01, tmax = 0. - i = 0 i = 1 i = 2 i = 3 i = - x = 0.00 x = 0.25 x = 0.50 x = 0.75 x = 1. - t = - n = 0 t = 0+ - n = 1 t = 0.01 - n = 2 t = 0.02 0 731.2 948.8 731.2 - n = 3 t = 0.03 0 649.0 879.2 649.0 - n = 4 t = 0.04 0 582.0 805.5 582.0 - n = 5 t = 0.05 0 524.6 734.0 524.6 - n = 6 t = 0.06 0 474.2 667.0 474.2 - n = 7 t = 0.07 0 429.2 605.3 429.2 - n = 8 t = 0.08 0 388.7 548.9 388.7 - n = 9 t = 0.09 0 352.1 497.7 352.1
- n = 10 t = 0.10 0 319.1 451.1 319.1
- n = 11 t = 0.11 0 289.2 408.9 289.2
- n = 12 t = 0.12 0 262.0 370.5 262.0
- n = 13 t = 0.13 0 237.5 335.8 237.5
- n = 14 t = 0.14 0 215.2 304.4 215.2
- n = 15 t = 0.15 0 195.0 275.8 195.0
- n = 16 t = 0.16 0 176.8 250.0 176.8
- n = 17 t = 0.17 0 160.2 226.5 160.2
- n = 18 t = 0.18 0 145.2 205.3 145.2
- n = 19 t = 0.19 0 131.6 186.1 131.6
- n = 20 t = 0.20 0 119.2 168.6 119.2
- Exact t = 0.20 0 125.1 176.9 125.1
- Error t = 0.20 0 5.8 8.2 5.8 - Table 3- - = 1, Δx = 0.25, Δt = 0.02, tmax = 0. Numerical solution of the Conduction Equation - i = 0 i = 1 i = 2 i = 3 i = - x = 0.00 x = 0.25 x = 0.50 x = 0.75 x = 1.
- n = 0 t = 0+
- n = 1 t = 0.02
- n = 2 t = 0.04 0 564.8 795.2 564.8
- n = 3 t = 0.06 0 457.9 647.7 457.9
- n = 4 t = 0.08 0 372.1 526.2 372.1
- n = 5 t = 0.10 0 302.3 427.6 302.3
- n = 6 t = 0.12 0 245.7 347.4 245.7
- n = 7 t = 0.14 0 199.6 282.3 199.6
- n = 8 t = 0.16 0 162.2 229.4 162.2
- n = 9 t = 0.18 0 131.8 186.4 131.
- n = 10 t = 0.20 0 107.1 151.4 107.1
- Exact t = 0.20 0 125.1 176.9 125.1
- Error t = 0.20 0 18.0 25.4 18.0
the calculation. This requirement for stability is separate from the requirement for accuracy. Accuracy depends on small sizes for Δt and Δx. We may want to use a larger step size for Δt to obtain a desired accuracy, however we cannot use a step size that violates the stability condition in equation [3-11] if use the explicit method.
The appendix to this section of notes contains a derivation for the truncation error for this approach. There are two equivalent equations for this truncation error. Both of these are written as an infinite series in equation [3-12].
2
2
2 1 2
2
1 2
k
n
i
k
k k
k
k
n
i
k
k k n k i
x
T
k f
k
k
t
x
T
k k
f
TE x [3-12]
We see that the truncation error depends in a complex way on the two finite-difference grid spacings, Δt and Δx. Choosing f = 1/6 eliminates the first (k = 2) term in the truncation error and gives a higher order truncation error.
The method just described is called the explicit method because all temperatures at the new time step can be written as an explicit function of temperatures at the old time step. Alternative, implicit approaches use more than one value of the temperature at the new time step in the finite- difference equation. The Crank-Nicholson method is the most popular implicit method. This method uses a finite-difference representation of the conduction equation at a time point midway between the two specified time grid lines. Although there is no computational point there, we can write this time derivative as follows, using a central-difference approximation. In this expression the step size is Δt/2, the distance from the location of the derivative to the location of each of the two points used in its computation.
2
1
2
2 2
1 2 2 1
1
[( ) ] [( ) ]
^
n
i
n i
n i
n i
n i
n
i x
T
O t
t
T T
O t
t
T T
t
T
[3-13]
Although this midpoint location gives us an easy way to compute the time derivative, there is no such easy way to compute the space derivative. We can use an average of the finite difference expression for the space derivative at the two time points, n and n+1. Such an average produces a truncation error that can be seen by adding the Taylor series for fi+1 and fi-1 from equations [2- 15] and [2-18] in the notes on basics.
4 ''''
2 ''' 1 1
3 '''
2 ' '' 1
3 '''
2 ' '' 1
h
f
h
f f f f
h
f
h
f f fh f
h
f
h
f f fh f
i i i i i
i i i i i
i i i i i
[3-14]
When we solve the last equation above for fi, we get the truncation error that arises when we use the arithmetic average of two neighboring points (fi+1 and fi-1) to approximate the value at the midpoint.
1 1 2
4 ''''
2 1 1 ''
Oh
h f f
f
h
f
f f
f i i i i i i i
[3-15]
Applying this result to the
2
1
2
2
n
x i
T
term gives the following result.
[( ) ]
1
2
2 2
2 2
1
2
2
O t
x
T
x
T
x
T
n
i
n
i
n
i
[3-16]
In this case the truncation error is proportional to [(Δt)^2 ] because we are taking an average of this derivative at two time points to represent the derivative at the midpoint. We can use the second- order equation for the second derivative for each term in this average. Doing this and combining equations [3-16] and [3-13] gives
[( ),( )] 0
2 2 2
1 1 2
1 1 1
1 1
1
2
1
2
2 2
1
^
O t x
x
T T T
x
T T T
t
T T
x
T
t
T
n i
n i
n i
n i
n i
n i
n i
n i
n
i
n
i
[3-17]
The finite difference that we obtain from equation [3-17] contains three temperatures at the new time point. We can rewrite this equation, introducing the factor f defined in equation [3-6] as follows.
i n in in T^ in Tin ^ f Tin
f
T
f
T f T
f
1 1
1 1
1 1
^ [3-18]
This equation describes a tridiagonal system of simultaneous linear algebraic equations. As we saw in the notes on basics, this set of equations is easily solved by the Thomas algorithm.
A final approach, known as the fully implicit approach, uses we use a backward difference for the time derivative and a central difference for the space derivatives at the point (xi,, tn+1). In this approach we write the two derivatives as follows.
[( )]
1 1 1
1 1
1
2
(^112)
O x
x
T T T
x
T
O t and
t
T T
t
T in in in
n
i
n i
n i
n
i
[3-19]
Substituting these finite difference expressions into the differential equation gives the following result.
[( ),( ) ] 0
2
1 1 1
1 1
(^11)
2
(^12)
O t x
x
T T T
t
T T
x
T
t
T in in in in in
n
i
n
i
[3-20]
Rearranging the finite-difference terms (ignoring the truncation error) and using our usual definition of f = αΔt/(Δx)^2 , gives the following result, after some algebra.
n i
n i
n i
n
fTi f T fT T
1 1
1 1
1 (^12 )^ [3-21]
Table 3- 4 Crank-Nicholson Solution of the Conduction Equation = 1, Δx = 0.01, Δt = 0.0005, tmax = 0. i = 0 i = 1 i = 2 i = 3 i = 4 x = 0 x = 0.01 x = 0.02 x = 0.03 x = 0. t = 0 1000 1000 1000 1000 1000 n = 0 t = 0+ 0 1000 1000 1000 1000 n = 1 t = 0.0005 0 - 73.35 423.96 690.85 834. n = 2 t = 0.001 0 352.75 305.27 440.73 599. n = 3 t = 0.0015 0 25.70 320.81 439.19 533. n = 4 t = 0.002 0 203.86 209.57 347.52 473. n = 5 t = 0.0025 0 56.79 252.91 334.12 422. n = 6 t = 0.003 0 141.46 177.47 298.20 397. n = 7 t = 0.0035 0 66.73 209.02 279.22 363. n = 8 t = 0.004 0 109.40 160.30 263.81 347. n = 9 t = 0.0045 0 68.71 179.63 245.6 8 324. n = 10 t = 0.005 0 90.79 148.20 237.92 311. n = 11 t = 0.0055 0 67.50 159.07 222.68 296. n = 12 t = 0.006 0 78.99 138.51 217.76 285. n = 13 t = 0.0065 0 65.08 144.07 205.56 273. n = 14 t = 0.007 0 70.94 130.31 201.68 264. n = 15 t = 0.0075 0 62.29 132.69 192.04 255. n = 16 t = 0.008 0 65.10 123.21 188.58 247. n = 17 t = 0.0085 0 59.50 123.75 180.95 241. n = 18 t = 0.009 0 60.65 117.00 177.71 234. n = 19 t = 0.0095 0 56.86 116.50 171.59 228. n = 20 t = 0.01 0 57.10 111 .53 168.52 222. n = 21 t = 0.0105 0 54.43 110.47 163.53 217. n = 22 t = 0.011 0 54.19 106.68 160.64 212. n = 23 t = 0.0115 0 52.22 105.35 156.49 208. n = 24 t = 0.012 0 51.73 102.36 153.78 203. n = 25 t = 0.0125 0 50.21 100.93 150.27 199. Exact t = 0.0125 0 50.43 100.66 150.48 199. Error t = 0.0125 0 0.216 0.272 0.212 0.
In addition to providing an accuracy comparison, these example calculations in Table 3-4 show that there is no inherent link between stability and physical reality. In these notes we derived the stability criterion for the explicit method by using arguments of physical reality. However, we see that the Crank-Nicholson method can produce intermediate results that are unrealistic physically, but the method still is stable. The fully implicit method does not produce such unrealistic results, but its accuracy is less than that of the Crank-Nicholson method.
The differences in the accuracy of the two methods increase as the calculation is run for a longer maximum time. Table 3-6 shows some overall results for the final time step when the calculations are done for a maximum time of 1. The errors in the fully implicit method, with its truncation error that is O[Δt, (Δx)^2 ] are significantly higher than those in the Crank-Nicholson calculation, which has a truncation error that is O[(Δt)^2 , (Δx)^2 ]. Both calculations took 2000 time steps with f = 5.0. Table 3-6 also shows the results of an explicit calculation that used 20, time steps to achieve the same maximum time. These time steps were the minimum allowed with the maximum value of f = 0.50. (If we wanted to use f = 1/6 to get the higher order truncation
error from the explicit method, we would need 60,000 time steps to reach the same maximum time.)
Table 3- 5 Fully-implicit Solution of the Conduction Equation = 1, Δx = 0.01, Δt = 0.0005, tmax = 0. i = 0 i = 1 i = 2 i = 3 i = 4 x = 0 x = 0.01 x = 0.02 x = 0.03 x = 0. t = 0 1000 1000 1000 1000 1000 n = 0 t = 0+ 0 1000 1000 1000 1000 n = 1 t = 0.0005 0 358.26 588.17 735.71 830. n = 2 t = 0.001 0 218.22 408.43 562.69 682. n = 3 t = 0.0015 0 166.26 322.13 460.74 578. n = 4 t = 0.002 0 139.05 272.65 396.35 507. n = 5 t = 0.0025 0 121.84 240.25 352.17 455. n = 6 t = 0.003 0 109.75 217.08 319.77 415. n = 7 t = 0.0035 0 100.65 199.49 294.81 385. n = 8 t = 0.004 0 93.50 185.57 274.85 360. n = 9 t = 0.0045 0 87.68 174.19 258.43 339. n = 10 t = 0.005 0 82.82 164.67 244.62 321. n = 11 t = 0.0055 0 78.69 156.56 232.81 306. n = 12 t = 0.006 0 75.13 149.54 222.55 293. n = 13 t = 0.0065 0 72.00 143.38 213.53 281. n = 14 t = 0.007 0 69.24 137.93 205.52 271. n = 15 t = 0.0075 0 66.77 133.05 198.35 262. n = 16 t = 0.008 0 64.55 128.66 191.88 253. n = 17 t = 0.0085 0 62.54 124.67 186.01 246. n = 18 t = 0.009 0 60.70 121.03 180.64 239. n = 19 t = 0.0095 0 59.02 117.70 175.71 232. n = 20 t = 0.01 0 57.47 114.62 171.16 226. n = 21 t = 0.0105 0 56.03 111.78 166.95 221. n = 22 t = 0.011 0 54.70 109.13 163.04 216. n = 23 t = 0.0115 0 53.46 106.67 159.38 211. n = 24 t = 0.012 0 52.30 104.36 155.96 206. n = 25 t = 0.0125 0 51.21 102.20 152.76 202. Exact t = 0.0125 0 50.43 100.66 150.48 199. Error t = 0.0125 0 0.779 1.542 2.273 2.
Table 3-6 also shows the execution time for calculations done on a 700 MHz Pentium III computer using a C++ program. The measured execution times are so short that they may not be meaningful. In general, solutions for this problem took about 0.5 to 0.6 microseconds per node per time step for the implicit and Crank Nicholson method and about 0.06 to 0. microseconds per node per time step for the explicit method. For example, a Crank Nicholson solution with Δx = 0.001 (1000 space nodes) and Δt = 0.00001 (100,100 time steps for a maximum time of 1) took 63.39 seconds. This is equivalent to (63.39 seconds) / (1000 space steps) / (100,000 time steps) ● (1,000,000 μsec/s) = 0.6449 microseconds per node per time step. To get a similar accuracy (rms error in temperature at the final time of 1x10-12) with the same Δx required a time step of 0.00000002 and an execution time of 4,110.4 seconds. This corresponds to a computation time of (4,110.4)(0.001)( 0.00000002) = 0.0822 microseconds per node per time step.
boundary conditions) we will have (N – 1)(M – 1) nodes at which we will have to solve an equation like [3-25]. We can illustrate a typical set of equations by considering a very small grid with N = 6 and M = 5. This sample grid is shown in Figure 1. In this case we have (6 – 1)(5 – 1) = 20 nodes
Figure 1. Example of small finite-difference grid.
We see that the form of equation [3-25] is to relate each node to its four nearest neighbors. The nodes related by the typical finite-difference equation are said to form a computational molecule as shown in figure 1. We can write this set of equations in a matrix form if we first choose an order for the equations. For example we can start in the lower left-hand corner and solve the equations for all the nodes in one row before moving on to the next row. In this case we start solving in the following order u 11 , u 21 , u 31 , u 41 , u 51 , u 12 , u 22 , u 32 , etc. When solving an equation near a boundary such as the equation for u 12 , the value of u 02 , is a known boundary condition. Since it is known, it can be moved to the right-hand side of the equation giving the following result.
(^1302) 2 12 22
2 11
2 u 2 u u u u [3-26]
Each of the four computational nodes located in a corner of the grid would have two known boundary values that would be placed on the right hand side of the equation. The equation for the computational molecule centered at node 11 would be written as follows.
(^1002) 2 12
2 11 21
2 2 u u u u u [3-27]
The entire set of matrix equations for the grid shown above is presented on the next page. For simplicity, the matrix equations use = 1 giving the general equation for x = y.
uij 1 ui 1 j 4 u (^) ij ui 1 j uij 1 0 [3-28]
i = 0 i = 1 i = 2 i = 3 i = 4 i = 5 i = 6 j = 5
j = 4
j = 3
j = 2
j = 1
j = 0
u 43
u 01 u 11
u 10
u 21
u 12
Boundary
nodes
u 51
u 44
u 53
u 42
u 33
Compu-
tational
Molecule
Figure 2. Matrix Representation of Second-Order Finite-Difference Equations for Laplace’s Equation with x = y
55 64
45
35
25
04 15
63
03
62
02
50 61
40
30
20
01 10
54
44
34
24
14
53
43
33
23
13
52
42
32
22
12
51
41
31
21
11
u u
u
u
u
u u
u
u
u
u
u u
u
u
u
u u
u
u
u
u
u
u
u
u
u
u
u
u
u
u
u
u
u
u
u
u
The coefficient matrix is a 20 by 20 matrix with 400 coefficients. Most of the coefficients are zero. All the nonzero coefficients are shown, but the large blank areas contain zeros that are not explicitly shown in this diagram. For Laplace’s equation the right hand side matrix contains only the boundary conditions where the value of u on a boundary has been specified.
Equation [3-32] can be simply coded as follows. As usual we assume that all the variables have been properly defined and initialized prior to this code. The variables Nx, Ny, and omega are used to denote the x and y grid nodes and the relaxation coefficient, ω. The unknown potentials, uij, are represented by the U[i][j] array.
for ( i = 2; i <= Nx; i++ ) for ( j = 2; j <= Ny; j++) U[i][j] = omega * (u[i][j-1] + u[i-1][j] + u[i+1][j] + u[i][j+1] ) / 4 +( 1 – omega ) * u[i][j];
To consider advanced iteration techniques we want to regard the set of equations to be solved, such as those in equation [3-28], as a matrix equation. An example of such an equation is shown in Figure 2. This appears confusing at first because we are converting equations on a two dimensional grid to a matrix notation where each grid node has an equation that provides one row of the matrix equation. However, even though we started with a one-dimensional grid, we still have a system of linear equations to solve and any such system can be represented by a matrix. To see this, examine the column of unknown potential values in the matrix in figure 2. Although these unknowns have two-dimensional subscripts, we could renumber these unknowns with a single subscript, giving the conventional matrix notation.
Once we recognize that the two-dimensional (or even-three-dimensional) finite-difference equations form the usual set of simultaneous linear equations that are represented by general matrix equation, Ax = b , we can use this matrix equation to analyze the process for solving this system by iteration. To do this we divide the A matrix into three different matrices.
A L D U sothat L D U x b [3-33]
Here L is a lower triangular matrix with all zeros on and above the principal diagonal. U is an upper triangular matrix with all zeros on and below the principal diagonal. D is a diagonal matrix. Although the division of A into L , D , and U in equation [3-33] is a general one, you should be able to see what L , D , and U would be for the A matrix in Figure 2.
The Jacobi iteration of equation [3-29] can be represented in a matrix form as follows;
Dx L Ux b x D L Ux D b ( 1 ) () ( 1 ) 1 ( ) 1 ( ) ( ) n n n n [3-34]
The second equation above shows that the matrix – D -1( L + U ) acts as an iteration matrix. Old values in x (n)^ are multiplied by this matrix to obtain new values in x (n+1). Matrix iteration processes like this one will converge if and only if the eigenvalue with the maximum absolute value in the iteration matrix is less than one.^4 This maximum eigenvalue is known as the spectral radius , ρ. Furthermore, the number of iterations, niter, required to reduce the error by a factor of 10-p^ is given by the following equation.
ln( )
ln( 10 )
ln( )
ln( 10 )
p n
p
iter^ [3-35]
For the Jacobi iteration process on Laplace’s equation with Dirichlet boundary conditions, the spectral radius, for large square grids with Nx = Ny = N can be approximated, for large N, by the following formula
(^4) A sufficient condition for convergence is that the solution matrix be diagonally dominant. This requires that
|aii| ≥ j i
aij. This is usually satisfied for finite difference equations.
for N N N N
Jacobi Jacobi large 2 2
ln( ) ln 1 2
(^1 )
2
2
2
2
2
(^)
[3-36]
When this equation for ρJacobi is substituted into equation [3-35] we find that the required iterations are approximately [2 p ln(10) N^2 ]/π^2. Thus, the number of iterations is proportional to the grid size squared. For the Gauss-Seidel iteration process, where we use the new values as soon as they are available, the iteration process can be represented by the following equation.
( ) ( ) [ ] ( 1 ) () ( 1 ) 1 ( ) L Dx Ux b x L D Ux b n n n n [3-37]
For this iteration matrix, -( L + D )-1 U , the spectral radius for the analysis of Laplace’s equation with Dirichlet boundary conditions, in the limit of large grid size, N, is given by the following equation.
for N N N N
GS^1 ln( GS ) ln^12 large
2
2
2
2
2
(^)
[3-38]
Here the number of iterations to reduce the error by a factor of 10-p^ is about half that required for the Jacobi iterations. The matrix form for simultaneous overrelaxation can be found, starting with a modification of equation [3-31]. In this modification, equation [3-31] produces an initial guess, x (n+1)*, that is subsequently corrected to obtain the final value for x (n+1). Thus we have the following start for the matrix representation of successive overrelaxation.
( ) [ ]
( 1 ) ( ) ( 1 )* ( ) ( 1 )* 1 () x x x x x L D Ux b n n n n n n where [3-39]
We can combine the two equations in [3-39] to eliminate x (n+1)*.
( 1 ) ( ) 1 ( ) ( ) ( ) [ ] n n n n x x L D Ux b x [3-40]
The initial value of x (n)^ in equation [3-40] can be premultiplied by the identity matrix, I , and the resulting equation may be solved to determine the iteration matrix.
x I L D U x L D b
( 1 )^1 ( ) 1 1 ( ) n n [3-41]
Here, the iteration matrix is (1 – ω) I – ( L + D )-1 U , and its eigenvalues will depend on the choice of the relaxation factor, ω. The optimal choice of the relaxation factor, for Laplace’s equation with Dirichlet boundary conditions and large grid size is given by the following equation.
2 1 1
2
Jacobi
optimum
[3-42]
Using equation [3-36] for ρJacobi gives the following result
Appendix – Truncation Error in the Numerical Solutions of the Diffusion Equation
To analyze the truncation error for the diffusion equation, we start by writing the truncation error as an infinite series. To derive this result we start with equation [2-28] from the notes on basics, which we used in the derivation of the finite-difference expression for the second-order derivative.
2 4
1 ^ 1
h
f
h
f i fi fi fi i [2-28]
We can rewrite this equation for a dependent variable from the diffusion equation, T. We use the usual two-dimensional grid notation to obtain a finite-difference equation for the second derivative in the diffusion equation. We will write all the terms in the infinite series for the truncation error.
2
2
2
2 2 2
2 1 1
k
x
x
x T
x
T
T T T
n k
k i
k
n k
i
n i
n i
n i
^ [3A-1]
Solving for the second derivative gives the usual second-order finite-difference equation with the truncation error expressed as an infinite series.
2
2
2 2
1 1 2
2
k
x
x
T
x
T T T
x
T k
k
n
i
k
n k i
n i
n i
n
i
^ [3A-2]
Similarly, we can start with the equation used to derive the first order forward difference derivative.
3 '''
2 ' ''
h
f
h
f i fi fih fi i [2-15]
Next, we rewrite this equation for potential as a series expansion in time, holding the spatial coordinate constant.
2
1
k
t
t
T
t
t
T
T T
n k
k i
k
n k
i
n i
n i
[3A-3]
We can solve this equation for the first derivative.
1
2
1
k
t
t
T
t
T T
t
T
n k
k i
k
n k i
n i
n
i
^
^ [3A-4]
We can substitute equation [3A-2] and [3A-4] into the diffusion equation to get the overall truncation error for replacing the differential equation by a finite-difference equation.
2 2
2
2
1 2
2
2
1 1
1 2
2
k
x
x
T
k
t
t
T
x
T T T
t
T T
x
T
t
T
k
k
n
i
k
n k k
k i
k
k
n i
n i
n i
n i
n i
n
i
n
i
[3A-5]
In equation [3A-5], the top row is the differential equation and the finite difference representation. The terms in the second row (except for the “=0”) are the truncation error. We can write this truncation error, denoted by the symbol, TE, at the given (i,n) grid point.
2 2
2
2
1 2
2 k
x
x
T
k
t
t
T
TE
k
k
n
i
k
n k k
k i
k
k n i
^
[3A-6]
By using the diffusion equation and the independence of the order of differentiation in mixed partial derivatives, we can derive the following result for higher order derivatives of the dependent variable, T.
4
4 2 2
2 2
2 2
2 2
2 2
2
x
T
x
T
t x
T
x x
T
t t
T
t t
T
[3A-7]
This leads to the following general result, which can be proved by induction.
k
k k k
k
x
T
t
T
2
2
[3A-8]
We can use this result to simplify equation [3A-6] for the truncation error.
2
2
(^12)
2
2 2 2
2
2 2 1 2
k
n
i
k
k k k k
n
i
k
k k k k n i
x
T
x k k
t
k
x
x
T
k
t
k
x
TE
[3A-9]
This equation shows that the final truncation error involves a relationship between the steps in
space and time and the infinite series for the truncation error contains the same 2
( x )
t
f
term
that appears in the finite-difference equation. We can write equation [3A-9] in terms of this ratio.
2
2
1 2 2 2
k
n
i
k
k k n k i
x
T
k
f
k
TE x [3A-10]
This equation hides the dependence of the truncation error on the time step in the f term. We could have done the second step in equation above analysis to obtain a result that uses only the time step and the f term. That would have given the following equation for the truncation error.
2
2
2 1
k
n
i
k
k k
k n i
x
T
k
k
k f
t
TE [3A-11]
