# E 2 u n1 1 2 u n ujn1 ujn j j 415 t x2 recall

This preview shows page 1. Sign up to view the full content.

This is the end of the preview. Sign up to access the rest of the document.

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 Crank-Nicolson method. We need a method of solving the initial-boundary 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 N-dimensional 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 space-e 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.

Ask a homework question - tutors are online