Let us now investigate problems having neumann

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: (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 higher-order 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 J-1 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 one-level 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.

Ask a homework question - tutors are online