Unformatted text preview: (4.2.1c) t>0 (4.2.1d) which has a Neumann condition at x = 1.
It is possible to use backward or forward di erences to discretize a Neumann boundary
condition at the right or left end of the domain, respectively however, it is generally
preferable to handle such boundary conditions by introducing a ctitious external grid
line and using higherorder centered di erences. For the problem at hand, construct a
uniform grid having cell spacing x = 1=J by t on 0 x 1, t > 0 and introduce the
ctitious grid line x = (J + 1) x as shown in Figure 4.2.1. With the extended problem
domain, we may approximate the derivative in the boundary condition (4.2.1d) using the
usual centered di erence formula
un ; un
(ux)n = J +1 x J ;1 + O( x2):
(4.2.2)
J
2
Neglecting the local discretization error in (4.2.2) and substituting the result into (4.2.1d),
we obtain an approximation of the Neumann condition at (j x n t) as
n
n
UJ +1 = UJ ;1 + 2 xgn: (4.2.3) Since centered di erences (4.2.2) were used, (4.2.3) is an O( x2 ) approximation of the
boundary condition (4.2.1d). This is consistent with the order of the centered nitedi erence approximations (2.1.9) that were used to discretize the partial di erential
equation (4.2.1a). 4.2. Neumann Boundary Conditions 13 t, n n 2
1
x, j
0 1 2 j J1 J J+1 Figure 4.2.1: Computational grid showing ctitious grid line at x = (J + 1) x to handle
a Neumann boundary condition.
The discrete problem on n t < t (n + 1) t is obtained by combining (4.2.3)
with the di erence approximation of the partial di erential equation, e.g., (4.1.7b), for
j = 1 2 : : : J . This gives a system of J + 1 equations for the J + 1 unknowns Ujn+1 ,
n+1
j = 1 2 : : : J + 1. It is often easy to eliminate the ctitious unknown UJ +1 from the
system. For example, suppose that we are discretizing (4.2.1a) using the weighted average
scheme (4.1.7b), then writing (4.1.7b) at j = J gives
+1
+1
;r UJn;1 + (1 + 2r )UJn+1 ; r UJn+1 = r(1 ; )UJn;1+
n
n
1 ; 2r(1 ; )]UJ + r(1 ; )UJ +1:
n+1
n
Now, eliminate UJ +1 and UJ +1 using (4.2.3) to get +1
;2r UJn;1 + (1 + 2r )UJn+1 = 2r(1 ; )UJn;1 + 1 ; 2r(1 ; )]UJn +2r(1 ; ) xgn + 2r xgn+1: (4.2.4) When (4.2.4) is combined with (4.1.7b) at the interior points j = 1 2 : : : J ; 1, we
obtain a system of J equations for the J unknowns Ujn+1 , j = 1 2 : : : J . This system
is tridiagonal and may be solved using the tridiagonal algorithm (4.1.14, 4.1.15). 14 Parabolic PDEs Example 4.2.1. Let us analyze the stability of the weighted average scheme (4.1.7b,
4.2.4) for solving the heat conduction problem (4.2.1) with f (t) = g(t) = 0 in order to
determine if a homogeneous Neumann boundary condition alters the stability characteristics of (4.1.7). We will use matrix methods since, in principle, the von Neumann method
ignores boundary conditions. (As with the Fourier series solution of partial di erential
equations, however, certain homogeneous boundary conditions may be treated without
di culty by periodic extension of the domain and initial data.)
Using (4.1.7b) and (4.2.4) with f (t) = g(t) = 0 yields the linear system I ; r C]Un+1 = I + r(1 ; )C]Un
2
6
Un = 6
6
4 U1n
Un
2 ...
n
UJ 3
7
7
7
5 2
;2 1
6 1 ;2 1
6
...
C=6
6
6
4
1 ;2 1
2 ;2 (4.2.5a)
3
7
7
7:
7
7
5 (4.2.5b) Equation (4.2.5a) may be written in the standard form Un+1 = L Un
with L = (I ; r C);1 I + r(1 ; )C]: We'll seek to bound the spectral radius of L . Thus, consider the eigenvalue problem
(L ; I)z = 0
or
or I + r(1 ; )C ; (I ; r C)]z = 0
r(1 ; + )C ; ( ; 1)I]z = 0: Write the above equation as (C ; I)z = 0 4.2. Neumann Boundary Conditions
where 15
1
= r(1 ; ;+ ) : Thus, is an eigenvalue of C and
r
= 1 +1 ; (1r; ) = 1 + 1 ; r r :
The eigenvalues and eigenvectors of many tridiagonal systems are known and, in particular, those of C are
2k
zjk = sin k2Jj
j = 1 2 ::: J
k = 2m ; 1
k = ;4 sin
4J
m = 0 1 : : : J ; 1:
(We'll indicate how to obtain these results in Chapter 9.)
Knowing k , we calculate k as
4r sin2 k =4J :
k =1;
1 + 4r sin2 k =4J
This is identical to the result (4.1.17) that we obtained for the ampli cation factor of
the homogeneous Dirichlet problem (4.1.7) using the von Neumann method thus, the
homogeneous Neumann condition has not altered the basic stability of the scheme.
The necessary stability condition that (L ) 1 is also su cient in this case since
L is similar to a symmetric matrix (cf. Mitchell and Gri ths 15] or Section 3.3. In
order to see this, let
3
2
1
7
61
7
6
7:
...
D=6
7
6
7
6
4
1p5
2
Then
2
3
;2 1
6 1 ;2 1
7
6
7
6
7:
...
^ = D;1CD = 6
C
p7
6
7
4
1 p2 2 5
;
2 ;2
is symmetric. Now let
^
L = D;1L D = D;1(I ; r C);1 I + r(1 ; )C]D 16 Parabolic PDEs or ^
L = D;1(I ; r C);1DD;1 I + r(1 ; )C]D
= D;1(I ; r C)D];1 D;1(I + r(1 ; )C)D]: ^
Using the similarity of C and C,
^
^
^
L = (I ; r C);1 I + r(1 ; )C]:
^
^
Taking the transpose of L and utilizing the symmetry of C, we nd
^
^
^
^
^
LT = (I + r(1 ; )C)T (I ; r C);T = (I + r(1 ; )C)(I ; r C);1 :
Now, ^
^
^
^
(I + r(1 ; )C)(I ; r C)];1(I + r(1 ; )C)(I ; r C) = I ^
^
however, I + r(1 ; )C and I ; r C commute, so
^
^
^
^
(I ; r C)(I + r(1 ; )C)];1 (I + r(1 ; )C)(I ; r C) = I
or ^
^
^
^
(I + r(1 ; )C);1(I ; r C);1(I + r(1 ; )C)(I ; r C) = I: Thus, ^
^
^
^
(I ; r C);1(I + r(1 ; )C) = (I + r(1 ; )C)(I ; r C);1 ^
^
and L = LT . 4.3 Multilevel Schemes and the Method of Lines
The onelevel schemes that we have been studying only use data at tn to obtain a solution
at tn+1 . Multilevel schemes employ solutions at tn and prior times to...
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