This preview shows page 1. Sign up to view the full content.
Unformatted text preview: pi = p(xi ), etc. Referring to this as the centraldi 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 pseudoPASCAL 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
 Spring '14
 JosephE.Flaherty

Click to edit the document details