This preview shows page 1. Sign up to view the full content.
Unformatted text preview: j
=
:
(4.1.5)
t
x2
Recall (Section 2.1) that the central di erence operator satis es Ujn = Ujn+1=2 ; Ujn;1=2 (4.1.6a) 2 U n = U n ; 2U n + U n :
j
j +1
j
j ;1 (4.1.6b) thus,
The parameter is a weighting factor selected between zero and unity. The computational
stencil for (4.1.5) when 0 < < 1 is shown in Figure 4.1.3. Several schemes are contained
within this weighted scheme. In fact, 4.1. Implicit Di erence Methods 5 if = 0, (4.1.5) is the explicit scheme (4.1.2),
if = 1, (4.1.5) is the backward Euler scheme (4.1.4), and
if = 1=2, (4.1.5) is an implicit scheme called the CrankNicolson method.
We need a method of solving the initialboundary value problem (4.1.1) when using
(4.1.5) with 6= 0. As usual, assume that the solution Ujn , j = 0 1 : : : J , is known and
that we want to determine Ujn+1 , j = 0 1 : : : J . Using the boundary conditions (4.1.1c)
we have U0n+1 = f n+1 n
UJ +1 = gn+1 (4.1.7a) where f n f (n t), etc. Using (4.1.6b) to evaluate the di erences in (4.1.5) and writing
the result at all interior mesh points j = 1 2 : : : J ; 1 gives
+1
+1
;r Ujn;1 + (1 + 2r )Ujn+1 ; r Ujn+1 = r(1 ; )Ujn;1+ 1 ; 2r(1 ; )]Ujn + r(1 ; )Ujn+1 j = 1 2 : : : J ; 1: (4.1.7b) n
If we regard U0n+1 and UJ +1 as being known quantities, then (4.1.7b) is a linear
tridiagonal system of J ; 1 equations for the J ; 1 unknowns Ujn+1 , j = 1 2 : : : J ; 1.
It is convenient to write (4.1.7b) in vector form, so let Un = U1n U2n : : : UJn;1]T (4.1.8) be the vector of unknowns at time level n. Then, write (4.1.7b) as I ; r C]Un+1 = I + r(1 ; )C]Un + rf n
where
2
;2 1
6 1 ;2 1
C=6
6
...
4
1 ;2 3
7
7
7
5 2 n+1
f + (1 ; )f n
6
0
6
...
n=6
f6
6
4
0
gn+1 + (1 ; )gn (4.1.9a)
3
7
7
7:
7
7
5 (4.1.9b) 6 Parabolic PDEs Tridiagonal systems, such as (4.1.9), arise frequently in numerical analysis and it
is important to have a fast algorithm for solving them. Fortunately, a special case of
Gaussian elimination su ces. Let us describe the method for a slightly more general
problem than (4.1.9) thus, consider a procedure for solving linear systems of the form AX = F (4.1.10a) where F and X are Ndimensional vectors and A is an N N tridiagonal matrix having
the form
2
3
a1 c1
6 b2 a2 c2
7
6
7
6
7:
...
A=6
(4.1.10b)
7
6
7
4
5
bN ;1 aN ;1 cN ;1
bN aN
In the present case, N = J ; 1, aj = 1 + 2r j = 1 2 ::: J ;1 (4.1.11a) bj = ;r j = 2 3 ::: J ;1 (4.1.11b) cj = ;r j = 1 2 ::: J ;2 (4.1.11c) and Fj = r(1 ; )Ujn;1 + 1 ; 2r(1 ; )]Ujn + r(1 ; )Ujn+1 + r
r j J ;1 gn+1 + (1 ; )gn] j1 f n+1 + (1 ; )f n]+ j = 1 2 : : : J ; 1: (4.1.11d) The quantity j k is the Kronecker delta which is unity when j = k and zero otherwise.
We assume that pivoting is not necessary (which is often the case when solving
parabolic problems by nite di erence or nite element methods) and factor A as A = LU: (4.1.12a) 4.1. Implicit Di erence Methods 7 Here, L and U are lower and upper bidiagonal matrices having the form
2
3
2
3
1
u1 v1
6 l2 1
7
6
7
u2 v2
6
7
6
7
6 l3 1
7
6
7:
...
L=6
U=6
7
7
6
6
7
... 7
4
5
4
uN ;1 vN ;1 5
lN 1
uN (4.1.12b) Once the coe cients lj , j = 2 3 : : : N , uj , j = 1 2 : : : N , and vj , j = 1 2 : : : N ; 1,
have been determined, the system (4.1.10a) may easily be solved by forward and backward
substitution. Thus, using (4.1.12a) in (4.1.10a) gives LUX = F: (4.1.13a) UX = Y (4.1.13b) LY = F: (4.1.13c) Let then, The system (4.1.13c) may be solved for Y by forward substitution and, once Y has been
determined, the system (4.1.13b) may be solved for X by backward substitution.
The coe cients of L and U are determined by a direct computation. Thus, substituting (4.1.10b) and (4.1.12b) into (4.1.12a), we nd u1 = a1
lj = bj =uj;1 uj = aj ; lj cj;1
vj = cj (4.1.14a) j = 2 3 ::: N j = 2 3 : : : N: (4.1.14b)
(4.1.14c) Having determined lj j = 2 3 : : : N , uj , j = 1 2 : : : N , and vj , j = 1 2 : : : N ; 1,
we can use (4.1.12b) in (4.1.13c, 4.1.13c) to nd Y1 = F1 Yj = Fj ; lj Yj;1 j = 2 3 ::: N (4.1.15a) 8 Parabolic PDEs Xj = (Yj ; vj Xj+1)=uj XN = yN =uN j = N ; 1 N ; 2 : : : 1: (4.1.15b) The procedure (4.1.14, 4.1.15) is called the tridiagonal algorithm. It can be implemented in a spacee cient manner by overwriting aj , bj , cj , and Fj with uj , lj , vj , and
xj (Figure 4.1.4). It requires about 2 divisions, 3 multiplications, and 3 additions or
subtractions for each j (other than j = 1 or N ). On some computers, division is much
slower than multiplication. In this case, we can trade 1 division for 2 multiplications by
storing 1=uj instead of uj . procedure tridi(N, a, b, c, F)
f Factorization g for j = 2 to N do
bj = bj =aj;1
aj = aj ; bj cj;1 end for f Forward and backward substitution g for j = 2 to N do
Fj = Fj ; bj Fj;1
end for FN = FN =aN
for j = N ; 1 downto 1 do
Fj := (Fj ; cj Fj+1)=aj end for Figure 4.1.4: Tridiagonal algorithm for solving (4.1.10).
Let us examine the stability of the weighted average scheme (4.1.7b) using the Von
Neumann method. As necessary, we assume that the initial conditions and solution are
periodic in x and write the solution as the discrete F...
View
Full
Document
This document was uploaded on 03/16/2014 for the course CSCI 6840 at Rensselaer Polytechnic Institute.
 Spring '14
 JosephE.Flaherty
 The Land

Click to edit the document details