313 in this case the approximation 6312a becomes yi1

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

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: pi = p(xi ), etc. Referring to this as the central-di erence approximation of the ODE we de ne the local discretization error as follows. De nition 6.3.1. Consider an ODE in the form Ly(x) = 0 and let Lhy be a di erence approximation of it, with L and Lh being di erential and di erence operators. The local discretization error or the local truncation error at x = xi is i = Lhy(xi): (6.3.15) Example 6.3.1. The di erential and di erence operators for the linear ODE (6.3.1, 6.3.13) satisfy Ly(x) = y00 + p(x)y0 + q(x)y ; r(x) = 0 and (x Lhy(xi) = y(xi+1) ; 2yh2 i) + y(xi;1) + p(xi) y(xi+1)2; y(i;1) + q(xi )y(xi) ; r(xi): h Using (6.3.7) and (6.3.8), we nd 2 h2 = y00(xi ) + 12 yiv ( i) + p(xi ) y0(xi ) + h y000( i)] + q(xi )y(xi) ; r(xi): i 6 Using the di erential equation 2 2 = h yiv ( i) + p(xi ) h y000 ( i): i 12 6 Thus, as we might have expected, the local discretization of the central di erence approximation of (6.3.1, 6.3.13) is O(h2). 15 The algebraic system (6.3.12b, 6.3.14) is linear for the linear BVP (6.3.1, 6.3.12b, 6.3.13) and may be solved by, e.g., Gaussian elimination. Towards this end, let us write (6.3.13) in the form biyi;1 + ai yi + ciyi+1 = h2 ri i = 1 2 ::: N ; 1 (6.3.16a) ci = 1 + h pi: 2 (6.3.16b) where bi = 1 ; h pi 2 ai = ;2 + h2 qi The boundary conditions (6.3.12b) may be used in conjunction with (6.3.16a) to create a system of dimension N + 1 or used to explicitly eliminate y0 and yN as unknowns from (6.3.16a). The latter approach is preferred for simple problems like this one. Thus, using (6.3.12b) with (6.3.16a) when i = 1, we nd a1 y1 + c1 y2 = h2r1 ; b1A: (6.3.17a) Similarly, using (6.3.12b) with (6.3.16a) when i = N ; 1 yields bN ;1 yN ;2 + aN ;1 yN ;1 = h2 rN ;1 ; cN ;1B: (6.3.17b) Grouping the N ; 1 equations, (6.3.17a), (6.3.17b), and (6.3.16a), i = 2 3 : : : N ; 2, yields Ay = f where (6.3.18a) 2 a1 c1 6 b2 a2 c2 6 A = 6 ... ... ... 6 6 4 bN ;2 aN ;2 cN ;2 bN ;1 aN ;1 2 6 y=6 6 4 3 y1 y2 7 7 ... 7 5 yN ;1 2 6 6 f =6 6 6 4 16 3 7 7 7 7 7 5 h2r1 ; b1 A h2r2 ... h2rN ;2 h2 rN ;1 ; cN ;1 B (6.3.18b) 3 7 7 7: 7 7 5 (6.3.18c) The linear algebraic problem (6.3.18) requires the solution of a tridiagonal system to determine the N ; 1 unknowns yi, i = 1 2 : : : N ; 1. The basic solution strategy is Gaussian elimination. Pivoting is frequently unnecessary. As seen from (6.3.16b), A will be diagonally dominant unless jp(x)j is large relative to jq(x)j or h is too small. Pivoting and other special treatment may be necessary in these exceptional situations. Let us proceed without pivoting and write A in the slightly more general form 2 3 a1 c1 6 b2 a2 c2 7 6 7 6 7 ... ... ... A=6 (6.3.19a) 7 6 7 4 bn;1 an;1 cn;1 5 bn an where, in our case, n = N ; 1. We factor A as A = LU (6.3.19b) where L is a lower triangular matrix and U is an upper triangular matrix. Having performed this factorization, we write (6.3.18a) as LUy = f let Uy = z (6.3.19c) Lz = f : (6.3.19d) and solve Since L is lower triangular, (6.3.19d) may be solved for z by forward substitution. Knowing z, we determine y by solving (6.3.19c) by backward substitution. All that remains is the determination of L and U. Let us hypothesize that they have the following bidiagonal forms 3 2 3 2 1 1 1 7 6 7 621 2 2 7 6 7 6 7: 6 7 6 ... ... 1 3 (6.3.20) U=6 L=6 7 7 7 6 6 ... ... 7 4 5 4 n;1 n;1 5 n n1 17 Using (6.3.20) with (6.3.19a,b) we nd = a1 (6.3.21a) = bi = i;1 (6.3.21b) 1 i i = ai ; i i = 2 3 ::: n i i;1 i = 1 2 : : : n ; 1: = ci (6.3.21c) (6.3.21d) Using (6.3.20a) and (6.3.19d) we nd the forward substitution step to be z1 = f1 zi = fi ; izi;1 (6.3.22a) i = 2 3 ::: n (6.3.22b) Similarly, using (6.3.20b) and (6.3.19c), the backward substitution step is yn = zn= yi = (zi ; iyi+1)= n i = n ; 1 n ; 2 : : : 1: i (6.3.23a) (6.3.23b) The solution procedure de ned by (6.3.19 - 6.3.23) is the basis of the famous tridiagonal algorithm. We state it as a pseudo-PASCAL algorithm in Figure 6.3.2. The version implemented in Figure 6.3.2 overwrites ai , bi , and ci with i, i, and i to reduce storage. By counting we see that the algorithm requires approximately 5n multiplications or divisions and 3n additions and subtractions. The work required to factor a full matrix by Gaussian elimination is approximately n3 =3. Thus, the ratio of the work to factor a tridiagonal matrix to that of a full matrix is approximately 15=n2. This is a significant savings even for small matrices and one should never use a Gaussian elimination procedure for full matrices to solve a tridiagonal system Example 6.3.2. Evidence from the Taylor's series expansion would suggest that the global error of the nite di erence solution of the BVP (6.3.1, 6.3.12b, 6.3.13) has an 18 Procedure tridi(n : integer var a b c f y : array 1::n] of real) begin f Factorization g for i = 2 to n do begin bi := bi =ai;1 ai := ai ; bici;1 end f Forward substitution g y1 := f1 fo...
View Full Document

Ask a homework question - tutors are online