New User Special Price Expires in

Let's log you in.

Sign in with Facebook


Don't have a StudySoup account? Create one here!


Create a StudySoup account

Be part of our community, it's free to join!

Sign up with Facebook


Create your account
By creating an account you agree to StudySoup's terms and conditions and privacy policy

Already have a StudySoup account? Login here


by: Addison Beer


Addison Beer
GPA 3.76


Almost Ready


These notes were just uploaded, and will be ready to view shortly.

Purchase these notes here, or revisit this page.

Either way, we'll remind you when they're ready :)

Preview These Notes for FREE

Get a free preview of these Notes, just enter your email below.

Unlock Preview
Unlock Preview

Preview these materials now for free

Why put in your email? Get access to more of this material and other relevant free materials for your school

View Preview

About this Document

Study Guide
50 ?




Popular in Course

Popular in Mathematics (M)

This 57 page Study Guide was uploaded by Addison Beer on Wednesday September 9, 2015. The Study Guide belongs to MATH 480 at University of Washington taught by Staff in Fall. Since its upload, it has received 28 views. For similar materials see /class/192107/math-480-university-of-washington in Mathematics (M) at University of Washington.

Similar to MATH 480 at UW

Popular in Mathematics (M)




Report this Material


What is Karma?


Karma is the currency of StudySoup.

You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 09/09/15
Chapter 11 Numerical Solution of TwoPoint Boundary Value Problems In the previous chapter we discussed initial value problems 2 ft7 t where 3tg is given For a single rst order equation like this the value of 3 at any one point determines its value for all time If 3 is a vector however so that we have a system of rst order equations or if we have a higher order equation such as uquota Fauau a 0 3 a 31 then more conditions are needed If the values of u and u are given at some point 30 E 01 then this ts the paradigm of an initial value problem but if instead say the values of u are given at a 0 and a 1 then this becomes a twopoint boundary value problem The restriction to the interval 01 is for convenience only the equation could be de ned over any interval a bl by making the change of variable a gt a b a If F is nonlinear in either or u the boundary value problem BVP is said to be nonlinear A second order linear two point BVP takes the form rau fa 0 3 a 31 111 u0 a u1 3 112 If is differentiable then the rst term can be differentiated to obtain puquot p39u but there may be good reasons to leave it in the form 111 111 An Application SteadyState Temperature Distribution Consider a thin rod as pictured below 0 X xAX 1 220 Suppose a heat source is applied to this rod and let ua 1 denote the temperature at position 3 time 1 By Newton s law of cooling the amount of heat owing from left to right across a in time At is HEE tAt where is the heat conductivity at position a in the rod thus heat ow is proportional to the gradient of temperature The net heat owing into the part of the rod between a and a An in time interval 1 to t At is therefore m Amiga Am 1m mag The net heat owing into this part of the rod in this time interval can also be expressed as an A a 7A1 a as c at where p is the density and 1A3 the mass of the section of rod 0 is the speci c heat the amount of heat needed to raise one unit of mass by one unit of temperature and 8tAt is the change in temperature over the given time period Equating these two expressions we nd 3 1A1 Am 81 1 fl 87 H83 quotmat When the system reaches a steadystate temperature distribution with at 0 this equation becomes 3 in 0 87 H83 This has the same form as equation 111 with H and E E E 0 where n which is now a function of a only 1 the st 1 1y st t 1 l If the ends of the rod are held at xed temperatures u0 a and u1 3 then we obtain the boundary conditions in 112 These are called Dirichlet boundary conditions If instead one or both endpoints are insulated so that no heat ows in or out at a 0 andor a 1 then the Dirichlet boundary conditions are replaced by homogeneous Neumann conditions u 0 0 andor u 1 0 Sometimes one has mixed or Robin boundary conditions 001mm cmu390 90 clgu1 cuu391 91 where egg cm cm 011 90 and 91 are given and taking the limit as As gt 0 gives 112 Finite Difference Methods Let us consider the special case where mm s 1 and qa s 0 in 111 uquota raua m 0 g a 31 113 um a W 3 114 221 Then it can be shown that if 2 0 then the problem 113 114 has a unique solution In general proving existence and uniqueness of solutions to boundary value problems is more dif cult than proving such results for initial value problems Given that the problem has a unique solution however one can then set about to approximate this solution numerically The rst step is to divide the interval 0 1 into small subintervals of width b 1n as pictured below i b0 x1 x2 xnillatn The mesh points or nodes are labeled 30 2 0 31 2 h 3 2 2h 37 2 nh 1 Next one approximates derivatives using nite difference quotients as described in Chapter 7 For example to approximate the second derivative uquot one might rst approximate u i h 2 using centered difference quotients Mitzi z ui1h uli z hui 1 and then combine these two approximations to obtain u ai h2 u ai h2 uai1 24331 h h Substituting this expression into 113 and letting n denote the value of the approximate solution at an we obtain the system of linear equations 1 ui1 2n 214 z 1 n 1 115 Combining this with the boundary conditions mg 2 a an 2 3 gives a system of n 1 equations in the n 1 unknown values ul un1 Moving all the known terms to the right hand side these equations can be written in matrix form as 2 h2ra1 1 U1 1 1 2 Warm 1 M2 1 2 h2rmn2 1 ung 1 2 h2TfEn1 un1 fa1 ah2 ak g 116 aw 2 I fltn 1 Notice that this system of linear equations is tridiayonal and can be solved with Gaussian elimina tion using only 0n operations See Section 524 222 1121 Accuracy As always we will be concerned with the accuracy of the approximate solution generated by such a nite difference method and especially in what happens in the limit as h gt 0 The global error or the co norm of the global error is de ned as 13113 M331 m where is the true solution and mi the approximate solution at node i The local discretization error or truncation error is the amount by which the true solution fails to satisfy the difference equations 1 rah E h fa 117 where is an exact solution of the differential equation It follows from 113 that 1 rah h uquota 118 The local truncation error can be analyzed using Taylor s theorem Assume that u E C39101 and expand h and h in Taylor series about 3 h h3 h391 h hu39a g quota 3721mm jaw 6 E lags h hi III hq IIII h hu a gu m 37a In r1 r1E h Substituting these expressions into 118 gives 1 h391 h391 I T3h hj EHWQ u39mm 0h2 Next we would like to relate the local truncation error to the global error This is usually the more dif cult part of the analysis We have from 115 and 117 lui1 2m 141 N331 an 0 i2ulti1 2mm alga 1 330M330 an an h Subtracting these two equations and letting ei E ui denote the error at Ode i we nd Vida 2 J 6141 330 WM i1n 1 119 Letting e E e1 en1T denote the vector of error Values these equations can be written in matrix form as Ac 2 T 223 where 39r E r1h quotran1 hT is the vector of truncation errors and 2 h2ra1 1 Ail 1110 1 2 h2ran1 is the same matrix as in 116 Assuming that A 1 exists the error is then given by eA 1r 1111 and we must study the behavior of A 1 as h gt 0 Assume rst that 2 y gt 0 for all a E 01 Recall that the in nity norm of a vector is the maximum absolute value of its components It follows from 119 that 2 TEih2 i i1 611 hgrmi h and hence that 2 WWmed S 2H6Hoc VHTHoc Since this holds for all i 1 n 1 it follows that 2 rhg lelloc S 2H6Hoc hgllTHx or Mac S HTHocV 0012 Thus the global error is of the same order as the local truncation error This same result holds when 0 but the proof is more dif cult Let us consider the case where E 0 There are various norms in which one can measure the error and now we will use something like the Lg norm 12 AWuW aw m where is the true solution and our approximate solution Unfortunately however the approximate solution is de ned only at mesh points This integral can be approximated by 71 1 12 wms g wm where denotes the usual Euclidean norm for vectors This is the composite trapezoid rule approximation to the above integral assuming that is a function that takes the values a and 3 at a 0 and a 1 respectively and u at a 2 mi 1 1 n 1 See Section 82 Note the factor As we re ne the mesh the number of mesh points 71 increases and so if we measured 224 the standard Euclidean norm of the error vector we would be summing the errors at more and more points When E 0 the matrix A in 1110 becomes A i quot quot 1112 h 1 12 which is a tridiayonal symmetric Toeplitz TST matrix It is tridiagonal because only the three center diagonals have nonzero entries it is symmetric because A 2 AT and it is Toeplitz because it is constant along any diagonal The main diagonal consists of 2h2 s the rst upper and lower diagonals of 1h2 s and all other diagonals of 0 s The eigenvalues and eigenvectors of such matrices are known Theorem 1121 The eigenvalues of the m by m TST matrix k7 A16 a 2bcos k 1m 1113 The corresponding orthonormal eigenvectors are k 2 V M7 211 lm sin k 1m 1114 Proof This theorem can be proved simply by checking that Gva Akva k 1 m and that vm 1 Since G is symmetric and has distinct eigenvalues its eigenvectors are orthogonal See Section 101 To derive the formulas however suppose that is an eigenvalue of G and v a corresponding eigenvector Letting 210 on 0 the equation Av v can be written in the form bw1 a w by 0 E 1 m 1115 This is a homogeneous linear difference equation of the sort discussed in Section 933 To nd its solutions we consider the characteristic polynomial b a 2 122 225 If C is a root of x then m CI is a solution of 1115 If the roots of x are denoted C1 and Cg then the general solution of 1115 is w chf ch if C1 and Cg are distinct or w chl cgEC if C1 Cg C where c1 and cg are determined by the conditions 210 on 0 The roots of x are a x a 4b2 2b A V A 24b2 1116 C1 2b C2 If these two roots are equal then the condition 210 0 implies that c1 0 and the condition 21m 0 implies that cgm 1Cm1 0 which implies that either cg 0 or C 0 and in either case 211 0 for all E Hence values of for which C1 Cg are not eigenvalues since there are no corresponding nonzero eigenvectors Assuming then that C1 and Cg are distinct the condition 210 0 implies that c1 cg 0 or cg 2 c1 Together with this the condition Um 0 implies that chi39l1 mg 2 0 which implies that either c1 0 in which case cg 0 and w 0 for all k or Cr 2 Cf The latter condition is equivalent to 2mk C2 C18X1gtlt7m1 1111 for some k 01 m We can discard the case k 0 since it corresponds to Cg C1 Substituting expressions 1116 for C1 and Cg into relation 1117 and multiplying each side by exp 2k7rm 1 we obtain the equation A a a 2 4b393exp Zk7r a a2 4b2exp m1 m1 Rearranging this becomes k A azsin m m 1 Squaring both sides gives 1 l 1 k7 1 k7 2 2 2 2 g 2 A a 4b cos 11 x a sin 11 l l 1 k7 2 4b 1 2 7 a os m 1 or Finally taking square roots gives ai2bcos k1m We need only take the plus sign since the minus sign gives the same set of values and this establishes 1113 226 Substituting 1113 into 1116 we readily nd that cos 2sin exp 21W m 1 in 1 in 1 k7 k7 zk7r 08 quotWm EXP W 1 M77 215k 2 01Cf 2012s1n 77711gt n H 9 and therefore Since this is a multiple of the vector va in 1114 it remains only to show that by choosing 01 so that va satis es 1114 we obtain a vector with norm one It is left as an exercise to show that m 2 M77 m1 Zsin 1 72 1 m and hence the vectors in 1114 have unit norm D Using Theorem 1121 we can write down the eigenvalues of the matrix A in 1112 1 k7 Ak 2 2cos k 1n 1 1118 Note that the eigenvalues of A are all positive A symmetric matrix with all positive eigenvalues is said to be positive de nite In order to estimate HeL2 using 1111 we must estimate HA lHL E max HA lWHL max lA11WWH2 EMA 1H2 2 med 2 Hmmwum Recall that the 2 norm of a symmetric matrix is the largest absolute value of its eigenvalues and that the eigenvalues of A 1 are the inverses of the eigenvalues of A Hence A 1 L2 is one over the smallest eigenvalue of A It follows from 1118 that the smallest eigenvalue of A is 1 1 2 2cos Expanding cos7rn cos7rh in a Taylor series about 0 we nd cos7rh 1 0h391 ltth 2 and substituting this into the expression for 1 gives 1 r2h2 00731 2 00 227 Hence A 1HL2 712 0h2 01 and it follows from 1111 that Helm S HA1HL2 HTHL2 01 CW 0072 VVe have thus established in the case where E 0 that the Lg norm of the global error is of the same order as that of the local truncation error ie 0h2 Note also from 1118 that the largest eigenValue of A is 2 2 WM Expanding cosn 17rn cos7r 7th in a Taylor series about 7139 gives 7711 T cos7r 7th 2 1 0h391 and it follows that 1 AH EM arth 0014 0012 Recall that the 2 norm condition number of a symmetric matrix is the ratio of its largest to smallest eigenValue in absolute Value Hence the condition number of A is Oh 2 This will be of importance later when we talk about solving linear systems iteratively 1122 More General Equations and Boundary Conditions Suppose in 111 is not identically 1 but is a function of 3 If gt 0 for all a E 01 and 2 0 for all a E 0 1 then it can again be shown that the boundary Value problem 1 M jg ru m 0 S a S 1 1119 110 a 111 3 1120 has a unique solution We can use the same approach to difference this more general equation Approximate dda by the centered difference formula 1 du z p dudmah2pdudaa h2 W h12u33 h UWDh a h2u33 W MW h 1 hjlwm W2 W h2u33 133 h2u33 h W h2u ML which is also second order accurate The difference equation at node i then becomes 1 h2 pai h2 pai h2ui1 pai h2ui1 fai h 228 and writing the difference equations at nodes i 1 n 1 in matrix form we obtain the linear system all b1 2L1 KMh 1 b1 a2 52 m fa2 hj bn 3 art 2 bn Q un l flmn Q I bn Q an l an fltn 1 bn lh2 where a pai h2 pai h2 h2rai i1n 1 b paih2 i01n 1 Note that the coef cient matrix for this problem is again symmetric and tridiagonal but it is no longer Toeplitz Still it can be shown that this matrix is also positive de nite Suppose now that in 111 is not identically 0 For simplicity let us return to the case where E 1 so that the differential equation becomes uquot qau39 fa 0 g a g 1 If we use a centered difference approximation to u z 1121 2h which is accurate to order hg then the difference equations become 1 u39 u39 Wm ui1 u11 maW mm mi 2 1n 1 This is the nonsymmetric tridiagonal system 01 b1 U1 f1 016 Cg a2 b2 Mg 122 E E 1122 Cn 2 tin 2 bn Q Uri 2 fltn 2 cn l an l un l fltn 1 n l where 2 1 q a 1 q a ai rmi bi ci E z1n 1 1123 For h suf ciently small the sub and super diagonal entries b and 0139 are negative and satisfy 2 1M 10139 229 If gt 0 for all i 1 n 1 then the diagonal entries a satisfy 2 a i Tquot 12139 139 Ci H h2 lt gt1b111 For this reason the matrix is said to be strictly row diagonally dominant De nition An m by m matrix A is row diagonally dominant if m km 2 2mm i1m 391 3 iii It is column diagonally dominant if m jam 2 2mm 391m 5 The matrix A is strictly row or column diagonally dominant if the corresponding inequality above is strict The importance of strict diagonal dominance is that it can be used to show that a matrix is nonsingular Clearly if the difference equations 1122 are to have a unique solution then the matrix there must be nonsingular It follows from Gerschgorin s Theorem Section 101 that a matrix that is strictly row or column diagonally dominant cannot have 0 as an eigenvalue and so must be nonsingular Hence for h suf ciently small the matrix in 1122 is nonsingular Sometimes the centered difference approximation 1121 to u is replaced by a one sided difference 1 z film u if lt 0 HQ ui1 1f 2 0 This is called upwind differencing It is only rst order accurate but it may have other advantages With this differencing 1122 holds with equations 1123 replaced by 2 a ai rmiqhill 1 11 39f 0 i if lt 0 b m h WWW h V 1 1 1 1 mango 1 IltT egt i qm20 z 7 Now the matrix is strictly row diagonally dominant when gt 0 for all i independent of h This is because bi W 2h2 while ai 252 As noted previously sometimes different boundary conditions are speci ed Let us return to the problem 113 114 Suppose the second Dirichlet boundary condition u1 3 is replaced by a Neumann condition u 1 3 Now the value an is no longer known from the boundary condition at a 1 but instead will be an additional unknown We will add an additional equation to approximate the boundary condition There are several ways to do this One idea is to approximate u 1 using a point outside the interval 1 uquot 1 z Then approximating the Neumann boundary condition by Mn 1 22m un1l 1 2 Mill 3 ElanH un 1l we can solve for an un1 un l 2W3 Substituting this into the approximate expression for u 1 and requiring that the differential equation 113 hold approximately at a 1 we obtain the equation hi2 This adds one more equation and one more unknown to the matrix equation in 116 1 2 2un1 2m mum mu 2 h2ra1 1 ul 1 1 2 h2ra2 1 ug 1 2 h2ran1 1 un1 2 2 h2ran an f aW f 1124 f mu 1 f mu 23M The matrix is no longer symmetric but it can be symmetrized Exercise 113 If in addition the rst Dirichlet condition u0 a is replaced by the Neumann condition u 0 a then ug becomes an unknown and a similar argument leads to the equation 1 hjpuu 2211 2ha T39fluuu Adding this equation and unknown to the system in 1124 we have 2 h2rau 2 W 1 1 2 h2TfE1 1 ul 1 2 h2ran1 1 un1 2 2 h2ran an f 330 20M f 331 E 1125 f mu 1 f 2 h If gt 0 for all i 0 1 n 1 n then this matrix is strictly row diagonally dominant and hence nonsingular On the other hand if 0 for all i then this matrix is singular because if one applies this matrix to the vector of all 1 s one obtains the zero vector This makes sense because the solution of the boundary value problem uquot f W0 a W1 is either nonexistent for example if f 0 but a y 3 or is nonunique if is any one solution then 7 is another solution for any constant 7 113 Finite Element Methods One can take a different approach to solving problems of the form 111w112 Instead of approx imating derivatives by nite difference quotients one might attempt to approximate the solution by a simple function such as a polynomial or piecewise polynomial While such a function might not be able to satisfy 111w112 exactly one could choose the coef cients so that these equations were approximately satis ed This is one of the ideas behind nite element methods Again let us restrict our attention for now to the BVP 113414 and let us assume that a 3 0 This problem can be written in the very general form Lu E uquot 0 g a g 1 1126 M0 u1 0 1127 In what sense should equation 1126 be approximately satis ed One idea is to force it to hold at certain points 31 mn1 This leads to what are called collocation methods Another idea is to choose the approximate solution say from the selected class of functions S in order to minimize the Lg norm of 22 f 12 Alma ue This leads to a least squares approximation Another idea is to note that if u satis es 1126 then for any function go the integral of the product of En with go will be equal to the integral of the product of f with go Let us de ne an inner product among functions de ned on 01 by 1 212051 uawada 232 Then if 21 satis es 1126 it also satis es 1 1 mums0 lt uquotltacrltmultasoltmdmA mamasum 1128 for all functions ga Equation 1128 is sometimes called the weak form of the differential equation One thing to note about the weak form 1128 is that we can integrate by parts in the integral on the left 1 uquots0 as u39s01 12212122wa 0 0 0 Hence 1128 can be replaced by 1 1 1 21 agaa110 2139aga ada0 1320 1129 With this formulation we can consider approximate solutions 21 that are differentiable only once instead of twice If 21 is restricted to come from a certain set S called the trial space then we probably will not be able to nd a function 21 E S that satis es 1128 or 1129 for all functions ga but perhaps we can nd a function 21 E S that satis es 1129 for all functions go in some set T called the test space The most common choice is T S Then taking S to be a space of piecewise polynomials this gives rise to a Galerkin nite element approximation More speci cally divide the interval 01 into subintervals of width b 121 with endpoints 30 2 0 31 2 h 3 nh 1 and let S be the space of continuous piecewise linear functions that satisfy the boundary conditions 1127 A function from S might look like that pictured below X0 X1 Xn We wish to nd the function 21 E S that satis es 1129 for all functions g0 E S To do this it will be helpful to have a basis for the space S a set of linearly independent functions in S with the property that every function in S can be written as a linear combination of these The set of hat functions pictured below forms a basis for S 91 p More precisely the set of continuous piecewise linear functions gai i 1 n 1 where g0 is equal to one at node i and zero at all other nodes ie 3 E mi 1321 megg1 i1n 1 1130 otherwise forms a basis for S To write an arbitrary continuous piecewise linear function go as a linear combination of these basis functions 21711 cigallm we simply take the coef cients 0 to be the values of the function go at each node Then the linear combination cm will have the same values as go at each of the nodes 31 mn1 as well as at the endpoints 30 and 37 where they are each zero and since a piecewise linear function is completely determined by its values at these nodes they must be the same function Example Suppose the interval 01 is divided into four equal subintervals and suppose is the continuous piecewise linear function de ned by 3 ifogmglzl 2m 14 imp13mg 12 in 74 2m ifl2gmgs4 1 3 if34gmgl The basis functions gal gag and gag for the space S of continuous piecewise linear functions that are zero at the endpoints of the interval have the value one at a 14 at a 12 and at a 34 respectively and they take the value zero at each of the other nodes Since 14 14 12 34 and 34 14 it follows that 14g01a 34gaga 14g03 To see that this is so note that 43 ingmgl4 4a 1 if14gmgl2 gala 2 43 if14gmgl2 gaga 3 43 if12gmg34 0 otherwise 0 otherwise 4m 2 if12gmg34 4 4a if34 3331 0 otherwise Hence 4ma raga4 mmltmgtltwmltwgtw4mlt3 E E Ei 3752 gram 4 4a1 ifs4931 Let 11s write our approximate solution as cjwjm The coef cients cl cn1 will be chosen so that Law 2 f gal for all i 1 n 1 The expression 21 gal should be understood to mean the integral after integrating by parts as in 1129 otherwise we would have trouble de ning 21 since 21 is discontinuous the second deriyative of 22 is zero in the interior of each subinterVal and in nite at the nodes Since any function go in S can be written as a linear combination of the gai s say 21711 dig and since the inner product is linear in each of its arguments this ensures that 21 f g0 2 22 f 21711 140 21711 di f gai 0 or equiyalently L ga f go More speci cally using 1129 we will choose cl cn1 so that 71 1 1 1 n1 1 a39ltmsoiltmlufu chso m s0 d0 rm Esme WW 11 11 lfagaiada i1n 1 0 This is a system of n 1 linear equations in the n 1 unknowns cl cn1 Note rst that the boundary term 21 agaa11 is zero since go in S implies 0 ga1 0 Note also that the summations can be pulled outside the integrals 71 1 1 1 1 2091 rim 20 13 1 1n 1 1131 11 The system 1131 can be written in matrix form as Ac 2 f 1132 where c E 01 cn1T is the vector of unknown coef cients f is the vector whose 1th entry is 01 13 i 1 n 1 and A is the n 1 by n 1 matrix whose ij entry is aij E f 13 f 13 1133 In other words the 1th entry of f is fgai and the ij entry of A is Lgaj gai after integrating by parts Returning to expression 1130 for g0 it can be seen that g0 satis es a E mi 1321 v a E miimiill 1134 otherwise It follows that aij is nonzero only for j i j 2 1 1 or j i 1 Hence the matrix A is tridiagonal It follows from 1133 that A is symmetric aij a for all ij 235 If in 1133 is a simple enough function then the integral there can be computed exactly Otherwise it can be approximated using a quadrature formula as described in Chapter 8 The same holds for the integrals in 1131 involving f Let us consider the simplest case in which E 0 Assume that the nodes are equally spaced 3 11 E h for all 2 Then 1 r5 1 2 n1 1 2 2 n 39 1 id 2 7 d f 7 7 f0 w a h w h h 1 1 1 1 1 NH am wiwi1 d3 h 3 d3 g 0 n5 If we divide each side of 1132 by h the linear system becomes 2 1 01 ltf 01gt 1 02 1 ltf1 02gt hi2 E f E 1 c Z Lian 2 1 2 cm mu 1 Note that the matrix on the left is the same one that arose in the centered nite difference scheme in 1112 Finite element methods often turn out to be almost the same as nite difference formulas Sometimes however what appear to be minor differences in the way boundary conditions are handled or in the way the right hand side vector is formed turn out to be important for accuracy Having derived the equations for the continuous piecewise linear nite element approximation for the problem 1126ww1127 it is easy to generalize to equation 111 Simply de ne the differential operator L appropriately and integrate by parts as before For equation 111 the i j entry of the coef cient matrix in 1132 would be if Ewy mz 1 1 1 0 dark0 dark0 d3 dgaj d3 1 1 1 01 1 f0 mmmmm da 1 1 1 0 mamaxmxazw da f0 qwxoxmxmm da 0 momma dm This matrix is still tridiagonal It is symmetric if E 0 236 1131 Accuracy One of the great advantages of nite element methods is that they sometimes permit easy proofs of the order of accuracy this is because the nite element approximation is optimal in a certain sense among all possible approximations from the trial space S When the operator L is selfadjoint ie 21717 2 2 217 for all 2 and w and positive de nite ie 2121 gt 0 for all 2 y 0 one might choose an approximate solution ii E S to minimize the energy norm of the error w 54 cm 54 wmm m mmmmm nxm Because L is self adjoint and positive de nite the expression in 1135 is a norm it often represents the energy in the system and sometimes the original problem is stated in this form Minimizing this expression leads to the Ritz nite element approximation Note that while we cannot compute L lf in 1135 since this is the solution that we are seeking we can compute the expression involving 22 on the right hand side and this expression differs only by a constant independent of 21 from the quantity on the left Hence it is equivalent to choose ii E S to minimize 2121 2f For the differential operator in 1119 this is f 321 a mama m 13 201 mmm 13 and after integrating by parts and using the fact that 210 2 211 2 0 this becomes 1 1 f0 mm W rltmltaltm 2gt da 20 mmmdm Expressing ii as djgaj we obtain a quadratic function in 11 dn1 2 1 71 1 2 71 1 1 71 1 dimidnifu pm Ema Arm deim ass 2A M Edisw in W quot39 11 11 11 The minimum of this function occurs where its derivative with respect to each of the variables 11 is zero Differentiating with respect to each 11 and setting these derivatives to zero we obtain the system of linear equations 1 71 1 I I 1 71 1 1 20211 ganja s01d2 f0 rm glam mam 2A fltacxailtmgtdmm LHWn L or equivalently n 1 i1r w7l 1 11 237 These are exactly the same equations that de ne the Galerkin nite element approximation This is no accident when E is self adjoint and positive de nite the Ritz and Galerkin approximations are the same Theorem 1131 IfL is selfadjoint and positive de nite then the Galerkin approximation is the same as the Ritz approximation for the problem Ln f Proof Suppose it minimizes 172 E v L lf 52 Elf over all functions 2 E S For any function 2 in this space and any number 5 we can write I21ev amp 1fev L1fev I2126 21 fv62 UU 1136 The only way that this can be greater than or equal to 172 for all e ie for 5 small enough in absolute value so that the 52 term is negligible and for 5 either positive or negative is for the coef cient of e to be zero ie 21 f 2 0 for all 2 E S Conversely if 21 f2 0 for all 2 E S then expression 1136 is always greater than or equal to I 21 so 21 minimizes I D This theorem is important because it means that of all functions in the trial space S the nite element approximation is best as far as minimizing the energy norm of the error Now one can derive results about the global error in energy norm just by determining how well can be approximated by a function from S eg by a continuous piecewise linear function say its piecewise linear interpolant Other trial spaces are possible too such as continuous piecewise quadratics or Hermite cubics or cubic splines Chapter 6 dealt with the errors in approximating a function by its piecewise polynomial interpolant so results from that chapter can be used to bound the error in the nite element approximation from a piecewise polynomial trial space The following theorem for example bounds the Lg norm of the error in the piecewise poly nomial interpolant of a function While the nite element approximation is not optimal in the Lg norm it can be shown that its order of accuracy in the Lg norm is the same as that of the optimal approximation from the space hence its order of accuracy is at least as good as that of the interpolant of u Theorem 1132 Let u be the k 1st degree piecewise polynomial interpolant of ii That is u matches u at each of the nodes mua1an1an and ifk gt 2 u also matches u at k 2 points equally spaced within each subinterval Then Hu 21er S ChklluklllL2 where h E mam 3111 mi and mm E f01va2 imp 3 238 For piecewise linear functions k 2 in the theorem and the theorem tells us that the Lg norm of the error in the piecewise linear interpolant is 0072 This is the order of the Lg norm of the error in the piecewise linear nite element approximation as well 239 Exercises 1 3 9quot 3 Consider the two point boundary value problem u 2mu 3221 32 u0 1 u1 0 a Let h 14 and explicitly write out the difference equations using centered differences for all derivatives b Repeat part a using the one sided approximation u tmi z h c Repeat part a for the boundary conditions u 0 u0 1 u 1 0 A rod of length 1 meter has a heat source applied to it and it eventually reaches a steady state where the temperature is not changing The conductivity of the rod is a function of position a and is given by 1 32 The left end of the rod is held at a constant temperature of 1 degree The right end of the rod is insulated so that no heat flows in or out from that end of the rod This problem is described by the boundary value problem fa 0 SL a Write down a set of difference equations for this problem Be sure to show how you do the differencing at the endpoints 1 Write a MATLAB code to solve the difference equations You can test your code on a problem where you know the solution by choosing a function that satis es the boundary conditions and determining what f must be in order for to solve the problem Try 1 Then 23352 23 1 c Try several different values for the mesh size h Based on your results what would you say is the order of accuracy of your method Show that the vectors va in 1114 are orthonormal ie show that m Iv lt7 1 wah 1 if j k otherwise Show that the matrix in 1124 is similar to a symmetric matrix Hintz Let D diag1 1 What is D lAD ifA is the matrix in 1124 240 5 Use the Galerkin nite element method with continuous piecewise linear basis functions to 05 solve the problem of exercise 2 but with homogeneous Dirichlet boundary conditions i 1x27 f7 0ltlt1 d3 d3 u0 0 u1 0 a Derive the matrix equation that you will need to solve for this problem 1 Write a MATLAB code to solve this set of equations You can test your code on a problem where you know the solution by choosing a function that satis es the boundary conditions and determining what f must be in order for to solve the problem Try 1 Then 2 2332 a 1 c Try several different values for the mesh size h Based on your results what would you say is the order of accuracy of the Galerkin method with continuous piecewise linear basis functions As noted in the text collocation methods choose an approximate solution from a space S so that the given differential equation holds at certain points in the domain Consider the two point boundary value problem uquota 352 0 3 a 1 a Consider the basis functions sinj7m j 12 3 and the collocation points 3 24 i 12 3 Suppose is to be approximated by a linear combination of these basis functions 21 of ja and suppose the coef cients are to be chosen so that the differential equation holds at the collocation points Write down the system of linear equations that you would need to solve to determine 0102 and cg b Repeat part a with ja 371 3 j 12 3 241 Chapter 12 Numerical Solution of Partial Differential Equations Previous chapters have dealt with differential equations in a single independent variable Chapter 9 dealt with time dependent problems while Chapter 11 dealt with problems having a single spatial variable In this chapter we will combine the two and deal with problems whose solutions depend on both time and space and or on more than one spatial variable These problems will be described by partial di erential equations A general second order linear partial differential equation in two independent variables a and 3 can be written in the form an 2m cum dufr euy in g 121 where the coef cients a b c d e f and g are given functions of a and 3 The equation is classi ed as one of three types based on the discriminant D E b2 ac If D lt 0 the equation is elliptic if D 0 it is parabolic and if D gt 0 it is hyperbolic Sometimes problems are of mixed type because the discriminant has one sign in one part of the domain and a different sign in another part Note that the lower order terms those involving u my and u are ignored in this classi cation The prototype of an elliptic equation is Poisson s equation um um g 122 Here b 0 and a 1 so that D E b2 ac 1l In the special case where g E 0 this is called Laplace 8 equation An elliptic equation with constant coef cients can be transformed to Poisson s equation by a linear change of variables Elliptic equations usually de ne stationary solutions that minimize energy The prototype of a parabolic equation is the heat equation also called the di usion equation u um 123 Here we have changed the variable name 3 to 1 because it usually represents time In the form 121 we have a 1 and b c 0 so that D 0 This equation was described in Section 111 242 where it was noted that uat might represent the temperature in a rod at position a and time t If there are two independent spatial variables a and 2 in addition to time t the heat equation becomes u um uw Parabolic equations usually describe smooth spreading ows such as the ow of heat through a rod in one spatial dimension or through say an electric blanket in two spatial dimensions A parabolic equation might also describe the diffusion of a liquid through a porous medium The prototype of a hyperbolic equation is the wave equation LL 2 mg 124 We have again used 1 instead of 3 since this variable usually represents time In the notation of 121 a 1 c 1 and b 0 so that D 1 Hyperbolic equations usually describe the motion of a disturbance preserving wave as for example in acoustics One classical example of a hyperbolic problem is a vibrating string such as a plucked guitar string Assume that the string is held xed at a 0 and a 1 Let ua 1 denote the displacement of the string at point a and time t and assume that the displacement is in the vertical direction only Let Tat denote the tension in the string ie the force that the part of the string to the right of a exerts on the part to the left of 3 at time 1 Assume that T is tangential to the string and has magnitude proportional to the local stretching factor x 1 as pictured below A Txt 0 x xh 1 If r is the tension in the equilibrium position along the horizontal axis then the magnitude of Ta t is quotrxl uat2 and the vertical component ofTat is quotrmAm 1 Now consider the part of the string between a and a h The vertical component of the force applied to this section of the string is TUlt h t TILAm 13 while the mass of this section of string is ph where p is the density and the acceleration is uLLat Hence by Newton s second law F 2 ma force equals mass times acceleration applied to this section of string 771 hit Tur t 1 h IAACE t Dividing each side by h and taking the limit as h gt 0 gives the equation Tumn 2 Wu and normalizing units so that p r we obtain the wave equation 124 121 Elliptic Equations We start with a study of time independent elliptic equations such as Poisson s equation in two space dimensions m yy 2 Lima 1 17 E Q 125 243 7 x 9 iyj 4 5 6 x y J Xi yj XH1 yj 1 2 X Y 3 1 r1 y0 x0 E hx E x1 Figure 121 Five point nite difference operator for Poisson s equation natural ordering of grid points 1 through 9 where Q is a region in the plane To uniquely de ne a solution we need boundary values say ua7 ga2 32 E 89 126 where 89 denotes the boundary of 9 1211 Finite Difference Methods Assume for simplicity that Q is the unit square 01 x 01 Divide the interval from a 0 to 1 into pieces of width hm 171 and divide the interval from 2 0 to 2 1 into pieces of width by 171 as in Figure 121 Let my denote the approximate value of the solution at grid point 3 394 E ihm jhy Approximate the second derivatives um and um by the nite difference quotients 2Uij 112414 iii 11 7 2Uij 2123941 TIM 1 2 hf h mrlmiz 3j z Hillmi yj z 127 There are 72 1ny 1 unknown values my i 1 n 1 j 1 ny 1 since the boundary values mm g03j unm g13j um gai 0 and um gai 1 are given If we substitute expressions 127 into equation 125 and require that the resulting difference equations hold at the interior nodes then we obtain 72 1ny 1 linear equations for the unknown values 2uv u39 x u39 2uv U39 U 11V l 1quot U 2 21 1quot 1fai3j z1n 1g1ny 1 2 hr 3 244 The nite difference operator on the left is sometimes called a vepoint operator because it couples my to its four neighbors 214151 2a um and 21314 It is represented schematically in Figure 121 In the special case where n my E n and hence hr by E 1n these equations become 41 11 lb 21 W 19 1 fai3j m 1 n 1 128 We can assemble these equations into a matrix equation but the exact form of the matrix will depend on the ordering of equations and unknowns Using the natural ordering equations and unknowns are ordered from left to right across each row of the grid with those in the bottom row numbered rst followed by those in the next row above etc as indicated in Figure 121 Thus point i j in the grid has index k i i j 1n 1 in the set of equations Using this ordering the matrix equation for 128 takes the form T I 111 f1 1h2u0 1 I T I 113 f2 129 I T I un3 fn2 I T u 1 111 1hun where T is the n 1 by n 1 tridiagonal matrix 4 1 T 1 1 1 4 I is the n 1 by n 1 identity matrix and um f 33139 Uj j0n f3 j1n 1 Uri 1 ak 1 31 This matrix is block tridiagonal because only the main diagonal blocks and the rst sub and super diagonal blocks are nonzero and it is symmetric It is also said to be black Toeplitz because the blocks along each block diagonal are the same Note however that while the individual blocks are Toeplitz the matrix as a whole is not because the rst elementwise sub and super diagonals contain 1 s except in positions that are multiples of n 1 where they contain 0 s The matrix in 129 is banded with half bandwidth n 1 It was shown in Section 524 that the work required to perform Gaussian elimination without pivoting on a banded matrix of size N with half bandwidth m is about 2ng operations During the process of Gaussian elimination zeros inside the band ll in with nonzeros which must be stored Hence the amount of storage required is about mN words For equation 129 N n 12 and m n 1 so the total 245 Figure 122 Differencing near the boundary of a general region 9 work required to solve the linear system using Gaussian elimination is about 2n 1 1 operations and n 13 words of storage are needed We will later discuss iterative methods for solving such linear systems which require less storage and sometimes less work than Gaussian elimination We also will see that the particular equations in 129 can be solved very rapidly using a fast Poisson solver As in the case of two point boundary value problems one can handle more general equations and more general boundary conditions using nite differences in two spatial dimensions One also can deal with more general regions 9 but this may become complicated Consider for example the region Q pictured in Figure 122 Suppose we wish to solve equation 125 on this domain with the same Dirichlet boundary conditions 126 Consider the difference equations at the point mugj in the gure Since 141304 lies outside the domain we cannot use it in our approximation of um at 31 Instead we would like to use the given boundary value which lies a distance say ah to the right of 3 along the grid line 3 394 Again we will use Taylor series to try and derive a second order accurate approximation to um at 31 Since all functions in our expansion will be evaluated at 2 394 we will simplify the notation by omitting the second argument First expanding 24331 about the point midj gives 2 3 um um mm 00 1210 246 Now expanding ah about the point mi 3 gives ah aha1a a 0011 1211 Adding or times equation 1210 to equation 1211 we nd a1ah2 owl Mb3 auai1 ah C l l ui Ui l 0011 and solving for um gives 1 39 1 39 h 1 1 39 1 h 1 mfgAmi 2 oauwz 1 11 a a u1 a Wm mi owl a1 ah2 3 Thus while one can approximate um at midj by using the grid values my and 214594 and the boundary value ah 394 2 orb3 the approximation is only rst order accurate the error is 1 ah3 mummy plus higher order terms in h In order to maintain a second order local truncation error one must use a difference formula that involves an additional point say 2 Expanding ua2 394 in a Taylor series about mi 3 and again dropping the second argument from the expression since it is always 394 we nd 2 ua2 2huai 7 21 001quot 1212 m 339 Let A B and C be coef cients to be determined and add A times equation 1210 and B times equation 1211 and C times equation 1212 Auai1 Buai ah Cua2 A B h A QB 2Cuai h2Aa2B4C hg Aa3B 8C 239 wnri 339 wmni In order for the left hand side divided by the coef cient of um on the right hand side to be a second order accurate approximation to um we must have AaB 2C0 Aa3B 8CO With a little algebra one nds the following second order accurate approximation to um u a luw39 1 Nb a2 14 cH l 11 a 1 aa1a2 ah A similar procedure must be carried out in order to nd a second order accurate approximation to 21W using values of u at mm942 mm941 mmj and mm4 h where 3h is the distance from 394 to the boundary along the grid line a Combining the two approximations one obtains the difference equation at node mi 247 1212 Finite Element Methods Finite element methods can be applied in two or more spatial dimensions as well Here we just brie y describe the procedure for solving a differential equation Ln f such as 125 with homogeneous Dirichlet boundary conditions ua 7 0 on 89 One can approximate ua 7 by a function ag from a space S called the trial space We will choose 22 so that 11 I if for all functions 23 in a space T called the test space Taking T S gives the Galerkin nite element approximation Here the inner product is de ned by 21217 EfQvazwazdmdy When evaluating 22 2 we will use the two dimensional analog of integration by parts namely Green s theorem ff mwm d3 d2 ff ayzay d3 dy 22123 17 Q Q 69 Here n denotes the derivative of 22 in the direction of the outward normal to the boundary 3939 Since we will be dealing with functions 1 that satisfy the homogeneous boundary condition 23 0 on 3939 however the boundary term will vanish As in the one dimensional case one must choose a trial space S and then nd a basis gal gaN of S One option is to take S to consist of continuous piecewise bilinear functions de ned on rectangles Suppose for example that Q is the unit square 01 x 01 One can divide the square into rectangles of width hr 171 in the m direction and of width by lny in the y direction The trial space would consist of functions that are continuous satisfy the boundary conditions and in each subrectangle have the form a bx 03 133 for some coef cients a b c and d A basis for this space can be formed by functions that are 1 at one node and 0 at all the others much like the hat functions in one dimension There is one such basis function for each interior node midj and a sample basis function is pictured below Each basis function is nonzero in only four rectangles The four parameters in its bilinear representation within a single rectangle are chosen to force it to have the desired values at each 248 of the four corners of the rectangle For example the basis function gaij1n21 associated with node Lj has the equation in mi 1334 gtlt zlj hzljl In mimii l gtlt zlj hzljl WHU mmx l m 1n 33433 gtlt 353041 in lii1l X BIN71 elsewhere Note that because the bilinear basis function is zero along the entire edges corresponding to a 351 3 2 an 3 2 3041 and 3 2 3941 it is continuous Another option is to use continuous piecewise linear functions de ned on triangles This has the advantage that a general region Q can usually be better approximated by a union of triangles than by a union of similar sized rectangles Again one can choose basis functions that are 1 at one node and 0 at all the others as pictured below Within each triangle a basis function has the form a bx c where the three parameters a b and c are chosen to give the correct values at the three nodes of the triangle A basis function associated with node i j is nonzero over all triangles that have node i j as a vertex As in the one dimensional case one writes the approximate solution 21 as a linear combination of the basis functions ag 2 ciga mn and then chooses the coef cients 01 cN so that N 11 The expression Lgaj gai is expressed using Green s theorem so that instead of involving second derivatives it involves only rst derivatives gain 94 gal and ga yl This is a system of N equations where N is the number of interior nodes in the grid in the N unknown coef cients which are actually the values of 22 at the interior nodes It has the form Ac 2 f where c 01 cNT f ltfs01gtm ltfa0vgt and Ala lt 0j2 0igt 122 Parabolic Equations We return to the time dependent heat equation in one spatial dimension L 2 up 1213 249 In order to uniquely specify a solution one must have initial conditions say at time t 0 as well as boundary conditions at the endpoints of the interval in 3 unless the equation holds in the in nite domain a E 00 For example if the equation holds for a E 0 1 and for t 2 0 then one might have the accompanying initial and boundary conditions ua0 qa 0 3 a 31 1214 u0t a u1t 3 t2 0 1215 where q is a given function of a and a and 3 are given constants Such a problem is called an initial boundary value problem IBVP 1221 Semidiscretization and the Method of Lines If we divide the spatial interval 0 1 into subintervals of width b 172 and let u1t un1t denote approximate solution values at each of the nodes 31 mn1 at time t then we can replace the spatial derivative um in 1213 by a nite difference quotient just as before 2u39t u39 tu39t unnit z A l 1 Substituting this approximation into equation 1213 we end up with a system of ordinary di er ential equations ODE s in time for the approximate solution values u1t un1t Q 2 4 t Re t iii 103 it h Letting ut denote the vector u1t un1t39l this system of ODE s can be written in the form u39 2 Au b where i1n 1 2 1 a 1 1 2 1 1 0 A2 bhj 1 2 1 0 1 2 3 Any of the methods discussed in Chapter 9 for solving the initial value problem for ODE s can then be applied to solve this system of ODE s with the initial value u0 qa1 qan1T The analysis of that section can be used to estimate the error in the numerically computed ap proximation to the solution of the ODE s while one also must estimate the difference between the true solution of the ODE s and that of the original PDE problem This approach of discretizing in space and then using an ODE solver in time is sometimes called the method of lines 1222 Discretization in Time Let us consider some speci c time discretizations Just as the spatial domain was divided into subintervals of width h let us divide the time domain into intervals of size At We will then approximate at using a nite difference in time and let n denote the approximate solution at 3 ih and tk kAt 250 Explicit Methods and Stability Using a forward di erence in time we approximate uat by uat At uatAt giving rise to the difference equations 211W ur on 2111 2112 At h Equivalently expressing the value at the new timestep in terms of the already computed values at the old timestep we have i1n 1 k01 1216 lcH lc At UE hj 2215 21931 21121 2 1 n 1 1217 where the boundary conditions 1215 provide the values ago 215 5 k01 1218 while the initial conditions give the values 21E qa i1n 1 1219 Thus formula 1217 gives a procedure for advancing from one time step to the next Starting with r 1 we can then compute the values 21 at time step two etc Using Taylor s theorem one can show that locally the method is second order accurate in space and rst order in time that is the local truncation error is 0h2 0At To see this substitute the true solution ua t into the difference equations 1216 and look at the difference between the left and right hand sides the values u we use formula 1217 with k 0 to obtain the values u 1 knowing those values uat At ua t 2ua t h t h t At h Expanding uat At about 3t gives ra t E ut At ut Atut t 0At2 and since 21 um this becomes aw At um t Ammo t 0am We have already seen that 2uat h t h t W 0W urrt and it follows that uat At mm At 2u t W 25quot W h t 0012 0mm 251 or subtracting ua t from each side and dividing by At t At t 2 t t t i 2 at at 747 uh 5 uh h 1220 am A It follows from 1220 that the local truncation error converges to 0 as h and At go to 0 but it does not follow that the global error the difference between the true and approximate solution converges to 0 To show that the approximate solution converges to the true solution on some time interval 071 is more dif cult and usually requires additional conditions on the relation between h and At The condition that rt gt 0 as h gt 0 and At gt 0 is a necessary condition for convergence and if this holds then the difference equations are said to be consistent In order to establish convergence of the approximate solution to the true solution one also must show that the difference method is stable There are a number of different notions of stability One de nition requires that if the true solution decays to zero over time then the approximate solution must do the same If we let a 3 0 in 1215 then since ua t represents the temperature in a rod whose ends are held xed at temperature 0 it follows from physical arguments that limtnx ua t 0 for all 3 Hence one should expect the same from the approximate solution values at the nodes 31 mn1 To see under what circumstances this holds we will write equations 1217 in matrix form and use matrix analysis to determine necessary and suf cient conditions on b and At to ensure that limknxui 0 for all i 1n 1 De ning 1 E Ath2 and still assuming that a 3 0 equations 1217 can be written in the form uk1 Aua where um E my u fll and 1 21 1 A quot 1221 1 1 1 21 It follows that ua Akuml and hence that ua converges to 0 as k gt 00 for all uml if and only if the matrix powers A C converge to 0 as k gt 00 This will be the case if and only if the spectral radius pA which is the maximum absolute value of an eigenvalue of A is less than 1 Fortunately we know the eigenvalues of the matrix A in 1221 According to Theorem 1121 they are E 1 2z2zcos l1n 1 For stability we need all of these eigenvalues to be less than 1 in magnitude ie 7 1lt1 2z 1 cosi lt1 1n 1 1222 n Subtracting 1 from each side and dividing by 2 this inequality becomes 1gt 11 1 cos gt 0 n 252 This holds for all E 1 n 1 provided 1 lt 1 2 since 1 cos 7rn is positive and is maximal when E n 1 and then its value is close to 2 if n is large Thus the necessary and suf cient condition for stability of the explicit method 1217 is At 1 7 A 1223 62 lt 2 l This is a severe restriction on the time step If h 01 for example this says that we must use a timestep At lt 1 20000 Thus to solve the problem for t between 0 and 1 would require more than 20000 time steps and would be very slow For this reason explicit time differencing methods are seldom used for solving the heat equation Example To see the effects of instability let us take say h 14 and At 14 so that 1 E Ath2 4 Starting with the initial solution 1 3 we set nag q14 316 as q12 14 and 2130 2 q34 316 Formula 1217 for advancing the solution in time becomes 21 421quot 4211 42112 2 1 2 3 k 01 239 21 and using this formula we nd 516 1916 18116 221116 u 14 u 34 ult3gt 594 u 4754 516 1916 18116 221116 After only 4 timesteps the co norm of 11 has grown from an initial value of 14 to 7754 3 194 and after 10 timesteps it is about 1e 9 The reason for this is clear The matrix 7 4 0 A 4 7 4 0 4 7 by which we are multiplying the approximate solution vector at each step has spectral radius 7 8cos37r4 m 1266 and so as we proceed in time this vector is being multiplied by approximately this factor at every step This is despite the fact that the true solution is decaying to 0 over time A plot of ul and U1 is shown in Figure 123 for the rst 4 time steps Another way to think of stability that leads to the same conclusion is to imagine that an error is introduced at some step which can be taken to be step 0 and to require that the error be damped rather than ampli ed over time Thus if instead of starting with the initial vector um one started with a slightly different vector 110 then the solution produced at step k would be 110 Ak m instead of um Akum The difference is um am 2 Aqua 130 and this difference approaches 0 as k gt 00 if and only if AC gt 0 as k gt 00 ie if and only if 1A lt 253 Figure 123 Unstable method h 14 At 14 Another approach to analyzing stability is sometimes called von Neumann stability analysis or the Fourier method We will now change notation and let the spatial Variable be denoted with subscript j since i will denote m For this analysis one often assumes an in nite spatial domain Assume that the solution a of the difference equations 1218 can be written in the form a gkemquot6 g cosjh isinjh 1224 where 9 represents an ampli cation factor since mil6 2 MI and g is a frequency in the Fourier transform of the solution Substituting this expression into 1217 gives glcleijh 1 mugceijh ngceij1ne Quayle from which it follows that g 1 21 1eih396 fl 1 2151 cosh The condition for stability is that be less than 1 for all frequencies g and this is the same as condition 1222 Hence we are again led to the conclusion that the method is stable if and only if 1223 holds Implicit Methods The nite difference method 1217 is called an explicit method because it gives an explicit formula for the solution at the new time step in terms of the solution at the old time step In an implicit 254 method one must solve a system of equations to determine the solution at the new time step This requires signi cantly more work per time step but it may be worth the extra expense because of better stability properties allowing larger time steps As an example suppose that the forward difference in time used to derive formula 1217 is replaced by a backward difference in time utat z uat uat AtAt The resulting set of difference equations is RikH RigQ 2 2ugk1 Rig1 Rig31gt At h i1n 1 k01 1225 Now to determine uk1 E u k1uki139r from um one must solve a system of linear equations Au k ultkgt b 1226 where 1 211 z a 0 Vquot A b 11 Z 1221 l 0 1 1 211 5 The local truncation error in this method is again 00 0At To analyze stability we again assume that a 3 0 so that b 0 in 1226 and ua A lua It follows that ua A kum and um will converge to 0 as k gt 00 for all um if and only if 1A1 lt 1 Once again Theorem 1121 tells us the eigenvalues of A E 121L 21LCOS 1n 1 and the eigenvalues of A 1 are the inverses of these 1 1n 1 1228 1 211 1 cos i i i Since 1 E Alth2 gt 0 and 1 cos 7Tn E 0 2 the denominator in 1228 is always greater than 1 Thus the eigenvalues of A 1 are positive and less than 1 for every value of 11 The method is therefore said to be unconditionally stable and it can be shown that the approximate solution converges to the true solution at the rate 00 0At independent of the relation between h and At Example When the implicit method 1225 is used with h 14 and At 14 the resulting approximate solution decays over time as it should Starting with 1 so that um q14q12q34T 31614316T the approximate solution vectors produced at time 255 0mm 0mm Figure 124 Stable method h 2 At 14 and h 18 At 116 Exact solution is plotted with dashed line steps 1 through 4 are 0 48 0103 0049 0015 um 0765 113 0230 ult3gt 0069 ult4gt 0021 0548 0103 0049 0015 The rst two components are plotted versus time in Figure 124 Also plotted with a dashed line is the exact solution at a 1 4 and a 1 2 While the method is stable it is not very accurate with h 14 and At 14 Also plotted in the gure are the approximate solution values using h 18 and At 1 16 Since the method is second order accurate in space and rst order in time the error in these values should be about 14 of that in the original approximation A potentially better implicit method is the CrankNicolson method which uses the average of the right hand sides in 1216 and 1225 H ugh ml 2 1 2u k uml u 2u k1 ugfjl air31 At 2 h b 1229 The average of these two right hand sides is an approximation to mfgAmi tC tk12 that is second order in both space and time that is one can show using Taylor s theorem that the error is 0h20At2 The left hand side of 1229 is a centered difference approximation to 21105 tci tk1 2 and so it is second order accurate that is the error is OAt2 Hence the local truncation error for the Crank Nicolson method is 00 OAt2 256 Moving terms involving the solution at the new time step to the left and those involving the solution at the old time step to the right formula 1229 can be written as 1 m1 1 M1 1 M1 1 1c 1 1c 1 1c 12 30 111 211 2 ui 1 l 111 211 2 ui 1 Once again we can analyze the stability of the method by using matrix analysis and Theorem 1121 or we can use the Fourier approach and substitute a Fourier component of the form 1224 into 1230 and determine whether such components grow or decay over time The two approaches again lead to the same conclusion Taking the latter approach and again letting j denote the spatial variable so that i can be used for V l we nd 1 glcleuh 1L9k16201h6 Em me 1ll gCeljh gllgklezmnme 620 1L or dividing by gkei f and solving for g 1 1 cosh 9 1111 cosh Since the term 11 cosh in the numerator and denominator is always positive the numerator is always smaller in absolute value than the denominator 1 z1 cosh lt 111 cosh 1z1 cosh Hence this method also is unconditionally stable Example Applying the Crank Nicolson method with the same values of h and At used in the previous examples yields the results in Figure 125 For h 2 At 14 the method is stable but not very accurate For h 18 At 1 16 the solution produced by the Crank Nicolson method can barely be distinguished from the exact solution plotted with a dashed line in the gure 257 mm Figure 125 Crank Nicolson method h 2 At 1 4 and h 1 8 At 116 Exact solution is plotted with dashed line 123 Separation of Variables We have discussed numerical solution of a number of model problems involving partial differential equations Actually some of these problems can be solved analytically using a technique called separation of variables This technique works only for problems in which the solution can be written as a product of two functions each of which involves only one of the independent variables In contrast the numerical techniques that we have described can be applied to general problems of the form 121 To solve the heat equation via separation of variables let us consider a problem with periodic boundary conditions utuw 033327139 tgt0 ua0 qa 0 3 a 3 27 u0t u27rt u0t 227139 t tgt 0 1231 This was a problem that Fourier studied It models the temperature in a circular rod of length 27139 To solve the IBVP 1231 via separation of variables assume that the solution ua t can be written in the form uat vawt for some functions 2 and M Then the differential equation becomes uaw t v awt 258 Assuming that y 0 and 217t y 0 we can divide each side by uawt to obtain w t 2 a 21713 Since the left hand side depends only on t and the right hand side depends only on 3 both sides must be equal to some constant say A Hence w t Awt tgt 0 1232 and 2 a Ava 0 3 a 3 27 1233 The general solution of equation 1232 is 2013 063 1234 Equation 1233 for 2 is an eigenvalue problem and must be coupled with the boundary conditions 210 2 21277 v390 v3927r 1235 in order for vawt to satisfy the boundary conditions in 1231 Equation 1233 for 2 always has 2 0 as a solution but its other solutions depend on the value of A We must consider three cases Case A lt 0 It can be shown in this case that the only solution to 1233 that satis es the boundary conditions 1235 is 2 0 Case ii A 0 In this case the general solution of 1233 is 2 c1 0251 The condition 210 2 21271 implies that c1 2 c1 cg27r or cg 0 Hence must be a constant and in this case the boundary condition 21 0 21 27139 holds as well Case iii A gt 0 In this case the general solution of 1233 is 2 c1 cosxrA cg sinxrA The boundary conditions 1235 imply that c1 2 c1 cosxrA 27 cg sinxrA 27 cg 2 c1 sinxrA 27139 cg cosxrA 27139 or equivalently limit gt i gt 3 gt This equation has a nonzero solution for c1 cg if and only if the matrix on the left hand side is singular that is if and only if its determinant is 0 1 cosxrA 271392 sin2xrA 277 0 In other words it must be the case that sinxrA 277 0 and 1 cosxrA 277 0 ie fo 1 2 Equivalently A mg m 1 2 In this case cosma and sinma are two linearly independent 259 solutions of the differential equation 1233 with the boundary conditions 1235 Substituting for A in 1234 and considering all of the cases iwwiii we nd that the solutions of the differential equation in 1231 that satisfy the periodic boundary conditions consist of all linear combinations of m2L 1 cosmae sinmaem2quot 22212 These are called fundamental modes of the differential equation with periodic boundary conditions To satisfy the initial condition ua0 qa we must expand in terms of these funda mental modes ie in a Fourier series 3C a0 2 am cosma bm sinma m1 A large class of 27r periodic functions q can be represented in this way In our problem q is de ned only on the interyal 0 27139 but assuming that q satis es the boundary conditions 10 2 q27r and q 0 q 27r it can be extended to be periodic throughout the real line Then assuming that the series converges appropriately the solution of the IBVP 1231 is x 2 21at a0 2 am cosma bm sinmaem I m1 One can check that this satis es the differential equation by differentiating the series term by term again assuming that the series converges in such a way that term by term differentiation is justi ed 3C 21 2 am cosma bm sin2m 7722em2l m1 3C 211M 2 Z am 2222 cosma bm 2222 sinmaem2quot 21 m1 Had the boundary conditions been different say homogeneous Dirichlet boundary conditions 210t 2122rt 0 instead of those in 1231 our analysis of problem 1233 would have been slightly different In this case the boundary conditions 1235 are replaced by 210 2 0 21221 0 1236 One again nds that in case i A lt 0 the only solution of 1233 that satis es the boundary conditions 1236 is 21 0 In case ii A 0 the general solution of1233 is again c cgm but now 210 2 0 01 0 and hence 21277 0 02 0 Thus 21 0 is the only solution in this case as well In case iii A gt 0 we again have the general solution 01 cosfo 02 sinxr 3 but the boundary condition 210 2 0 implies that 01 0 and the boundary condition 21277 0 then implies that sinxrA 277 0 or fo 12 Thus the fundamental modes for this problem are sinma m 12 260 Again to satisfy the initial condition one must expand in terms of these fundamental modes ie in a Fourier sine series This is the appropriate series to use for functions that satisfy the boundary conditions 10 2 q27r 0 Assuming that the series converges in the right way the solution to the IBVP is then 3C ua t Z bm sinmaem2l 1237 m1 The method of separation of variables works for a given problem provided the set of fundamental modes for the differential equation with the given boundary conditions is rich enough that the initial function can be represented as a linear combination of these fundamental modes Otherwise the solution ua 1 cannot be written as a product vawt and a different approach must be used to solve or approximately solve the IBVP We argued previously based on physical grounds that the solution to the heat equation in a rod with the ends held xed at temperature 0 would decay to 0 over time Now this can be seen analytically through expression 1237 each of whose terms approaches 0 as t gt 00 for all 3 1231 Separation of Variables for Difference Equations Just as separation of variables and Fourier series were used to solve the IBVP 1231 they can be used to solve the difference equations 1217ww1219 assuming that a 3 0 Assume that a can be written in the form 6 u impC Making this substitution in 1216 we have 39 r 2 r r 39 U17UIC1At 212016 2 212016 21122016 21 120k 1PH n 1 k 01 Collecting terms involving k on the left and those involving i on the right and assuming that 2 y 0 and wk 75 0 gives wk wk 2m 2141 251 NUg 71239 i where 1 Ath2 Since the left hand side is independent of i and the right hand side is inde pendent of k both sides must be equal to some constant say thus wk1 1 zwk k 01 1238 221 21139 m1 2 Am i1n 1 210 2 0 1239 261 Equation 1239 represents the eigenvalue problem for the TST matrix with 2 s on its main diagonal and 1 s on its sub and super diagonals From Theorem 1121 these eigenvalues are 2 2cos 7rn E 1 n 1 and the corresponding orthonormal eigenvectors are vm x2n sin 7rnsin2 7rnsinn 1e7rnquot 1240 Thus for each the vector v with 2 27 sini 7rn i 1 n 1 is a solution of 1239 The solution of 1238 is wk 2 1 Amkum and so the solutions of the difference equation 1216 with homogeneous Dirichlet boundary conditions are linear combinations of a vimC 1 21 211cos 7rnk2n sini 7rn l 1 n 1 To satisfy the initial conditions nip 2 qai i 1 n 1 one must expand the initial vector in terms of the eigenvectors in 1240 that is one must nd coef cients a1 E 1 n 1 such that n l Ea 27 sini 7rn gm i1n 1 1 Because the eigenvectors span fin 1 this expansion is always possible and because the eigenvectors are orthonormal it turns out that n l a1 x2anaj sinj 7rn E 1 n 1 11 Thus the solution of 1217ww1219 with a 3 0 is quot 1 n l my 2 211 21 21 cos 7rnk2n sinj 7rn sini 7rn 9 We can use this expression to answer questions about stability of the difference method as was done before The approximate solution gt will approach 0 as k gt 00 if and only if 1 21 21cos 7Tn is less than one for all E 1 n 1 We have seen previously that a suf cient condition for this is 1 g 12 124 Hyperbolic Equations 1241 Characteristics We begin our study of hyperbolic equations by considering the one way wave equation u an 2 0 a E coco tgt 0 ua0 1241 Here we consider an in nite spatial domain so no boundary conditions are present To uniquely specify the solution one still needs an initial condition ua0 qa so this is called an initial value problem IVP This problem can be solved analytically The solution is uat qa at 1242 262 as is easily checked by differentiating uat aq a at and uat q a at hence L 0 0 The solution at any time t gt 0 is a copy of the original function shifted to the right if a is positive and to the left if a is negative as pictured below F0 gt0 0 x 0 xat Thus the solution at 3t depends only on the value of q at a at The lines in the 3 1 plane on which a at is constant are called characteristics The solution is a wave that propagates with speed a without change of shape Note also that while the differential equation seems to make sense only if u is differentiable in 3 the solution 1242 is de ned even if q is not differentiable In general discontinuous solutions are allowed for hyperbolic problems An example is a shock wave which can occur in nonlinear hyperbolic equations Consider the more general equation uat auat buat fat a E coco tgt 0 1243 Based on the preceding observations let us change variables from 3 t to E r where Tt m at m ar De ne 215 7 uat Then 821 393th 8u83 3937 81187 8387 2 L awn since 3937 1 and 3937 a Hence equation 1243 becomes 821 E T 1727 This is an ordinary di erential equation ODE in r and the solution is T 2m r meM f m germ 9 ds 0 Returning to the original variables L uat qa awe bquot fa at s ska W s ds 0 Note that uat depends only on values of q and f at points 333 such that a af 2 a at ie only on values of q and f along the characteristic through 3 t This method of solving along characteristics can be extended to nonlinear equations of the form at an 2 fa t 263 1242 Systems of Hyperbolic Equations De nition A system of the form ut Aux Bu 2 fa t 1244 is hyperbolic if A is diagonalizable with real eigenvalues Example As an example of a hyperbolic system suppose we start with the wave equation u a um where a is a given scalar and we introduce the functions 2 2 an and w 2 up Then we obtain the system of equations 2 0 a 2 ina 0 217 1245 The eigenvalues of the matrix in 1245 are the solutions of X2 a 2 0 ie ia These are real and the matrix is diagonalizable The eigenvalues of A are the characteristic speeds of the system If A SAS I where A diag1 Am then under the change of variables 1 S lu equation 1244 becomes upon multiplying by S 1 on the left at A x 5 133 slf When B 0 these equations decouple 21 Mai S 1f i1m Example The eigenvectors of the matrix in 1245 are 227 corresponding to a and Mg2 x 2T corresponding to a Hence 5 22 f22 1 f22 f22 m2 2 m2 2 De ne 211 2v 2w and 212 52 u Then 221 221 so that 1mt My at and 212 a 2 so that 2mt 2 623 at Thus the solution is a combination of a wave moving to the left with speed a and another wave moving to the right with speed a 264 1243 Boundary Conditions Hyperbolic PDE s may be de ned on a nite interval rather than the entire real line and then boundary conditions must be speci ed This gives rise to an initial boundary value problem but one must be sure that the correct types of boundary conditions are given Consider again the one way wave equation m aufn 0 0 3 a 1 tgt 0 1246 de ned now on the nite interval 01 If a gt 0 the characteristics in this region propagate to the right while if a lt 0 they propagate to the left as illustrated below agt0 alt0 t t 0 x 1 0 x 1 It follows that in addition to the initial data at time t 0 the solution must be given on the left boundary of the region 0 when a is positive and on the right boundary 1 when a is negative If we specify initial data ua0 and boundary data u0t 2 gt for a gt 0 then the solution is qa at ifm atgt0 quot gt aa ifm atlt0 For a lt 0 the roles of the two boundaries are reversed 1244 Finite Difference Methods One can imagine any number of ways of differencing the equations 1246 Letting n denote the approximate solution at spatial point 394 2 jh and time 1116 kAt one could use a forward difference in time and a forward difference in space to obtain the equations MUM1 uEc uEc ulc J J a J 1 J At Alternatively one could try a forward difference in time and a backward difference in space 4 0 j1n 1k01 1247 M1 1c uj uj At af j1n 1k01 1248 A third alternative is a forward difference in time and a centered difference in space M1 6 k k u u u 1 u1 I W a gf 20 j1n 1k01 1249 265 The leapfrog scheme uses centered differences in time and space 16 k u u a H1 1 m1 lc l 9 9 J 2m 2h j1n 1k01 1250 The leap frog method is called a multistep method because it uses solutions from two previous time steps steps k and k 1 to compute the solution at the new timestep step k 1 In order to start the method one must use a one step method to approximate the solution at t1 given the initial data at time 110 0 Once the approximate solution is known at two consecutive time steps formula 1250 can then be used to advance to future time steps Still another possibility is the LamFriedrich scheme 906 l At 2h 0 j1n 1 k01 1251 The local truncation error LTE for these different methods can be analyzed just as was done for parabolic problems Using Taylor s theorem to determine the amount by which the true solution fails to satisfy the difference equations one nds that the LTE for the forward time forward space method 1247 is 0At 0h and the same holds for the forward time backward space method 1248 The LTE for the forward time centered space method 1249 is 0At 0h2 since the centered difference approximation to 2 is second order accurate The leap frog scheme 1250 has LTE that is second order in both space and time 0At2 001 The Lax Friedrichs scheme 1251 has LTE 0At 0072 0h2At with the last term arising because uaj1tk uaj1tc is an 0h2 approximation to uajtc but this error is divided by At in formula 1251 As we have seen before however local truncation error is only part of the story In order for a difference method to generate a reasonable approximation to the solution of the differential equation it also must be stable Stability and the CFL Condition The Lax Richtmyer equivalence theorem states that a difference method for solving a well posed partial differential equation is convergent if and only if it is consistent meaning that the local truncation error approaches 0 as h and At go to 0 and stable Without giving a formal de nition of stability we note that since the true solution to equation 1241 does not grow over time the approximate solution generated by a nite difference scheme should not grow either Consider for example the forward time forward space difference scheme 1247 If we de ne E Ath then 1247 can be written as k 1 k k k u lu au1 u 266 Looking at the sum of squares of the components on each side we nd 2 uk12 2 1 auk Mugif I 2 2 Z 1 my up 21 ammgmygzl W mg31 s 11aA1aAQZukgt2 9 A suf cient condition for stability is therefore 1a a g 1 or equivalently 1 g a g 0 It turns out that this is also a necessary condition for stability Since E Ath is positive the method cannot be stable if a gt 0 For a lt 0 the requirement is that g 1a or At 3 l 1252 la One can see in another way that condition 1252 is necessary for the difference method 1247 to be able to generate a reasonable approximation to the solution of the differential equation 1241 Consider the approximation my to the solution at time 1116 kAt and spatial point 394 2 jh We know that the true solution is uaj 116 2 qaj atC Suppose a is positive Now a is determined from 21906 1 and These values are determined from ask 2 ask 12 and Continuing to trace the values back over time we obtain a triangular structure of dependence as illustrated below x xh kh J J x39 Thus It depends only on values of the initial solution q between a 394 and a 394 kh Yet if a is positive then the true solution at 394 116 is determined by the value of q at a point to the left of 394 namely qaj akAt Hence despite the fact that the difference scheme is consistent it cannot generate a reasonable approximation to the true solution For a negative we still need the domain of dependence of my to include that of the true solution that is we need my kh to be greater than or equal to 394 akAt which is equivalent to condition 1252 267 The idea that the domain of dependence of the approximate solution must include that of the true solution is called the CFL ComantFriedrichsLewy condition More precisely we can state the following Theorem 1241 Consider a consistent di erence scheme of the form lc1 C k 6 l quill my yu l for equation 1241 A necessary condition for stability is the CFL condition Atgl 1a Proof The true solution uajtC is qaj atc On the other hand the approximate solution a is determined by iii 11 uk1 and 2641 these in turn are determined by 4132 etc Continuing it follows that my depends only on values of q between 394 kh and 394 kh Hence in order for the method to be convergent it is necessary that my akAt E kh 394 kh This means that kh 2 akAt or At 3 ha By the Lax Richtmyer theorem a consistent method is convergent if and only if it is stable Hence the condition At 3 ha is necessary for stability 1 lc 239 In addition to using the CFL condition to determine necessary conditions for stability the stability of difference methods for hyperbolic problems can be analyzed in a way similar to what was done for parabolic problems namely by assuming that the approximate solution has the form gkei f for some ampli cation factor 9 and some frequency g The condition for stability is g g 1 for all frequencies Using this approach let us consider the forward time centered space method 1249 Substituting g ceijh6 for my and dividing each side by gkeijh6 gives 1 H16 ih L a 0 At 2h or with Ath g 1 9er e M 1 ai sinh Since g2 1 a sin hf 2 gt 1 this method is always unstable Let us again consider the two way wave equation u alga 0 g a g 1 t gt 0 1253 It was shown in the example of Subsection 1242 that this second order equation could be reduced to a system of two rst order hyperbolic equations One way to solve this problem numerically would be to carry out this reduction and then apply an appropriate difference method to the system of rst order equations Another approach is to difference the equation as it stands In 268 addition to the differential equation one needs initial conditions and if the spatial domain is nite appropriate boundary conditions in order to de ne a unique solution Since the equation involves a second derivative in time we will need two initial conditions such as ua0 qa ua0 pa 0 3 a 31 If we take homogeneous Dirichlet boundary conditions u0t u1t 0 then for stability we again require that the approximate solution not grow over time Suppose a centered time centered space method is used for equation 1253 6 16 1 161 16 k k 2uj 2194 2194 a2 2uj uj1uj1 At h 1254 This is called a multistep method because the approximate solution at time step k 1 depends on that at step If as well as that at step k 1 To start the solution process we need to know um which can be obtained from the initial condition nip 2 qaj but we also need um This can be obtained from the other initial condition ua 0 Approximating uaj 0 by a u 0At we nd a 2 2190 Atpaj With the values of 11 known at timesteps 0 and 1 equation 1254 can then be used to advance the solution to timestep 2 etc Since centered differences are used in 1254 for both nu and um the local truncation error is 0At2 0h2 One can study the stability of the method by assuming that my has the form gkei and determining conditions under which g 1 Making this substitution in 1254 and dividing by g cei f gives 2 g 9 1 2 2 6mg 6 ih396 my a h or with Ath 2 g 9 1 a22 2 2cos hf After some algebra it can be shown that g 1 if and only if a g 1 Hence we are again led to the stability condition 1252 125 Fast Methods for Poisson s Equation We have seen that the solution of elliptic equations and parabolic equations often involves solving large systems of linear equations For Poisson s equation this must be done once while for the heat equation if an implicit time differencing scheme is used then a linear system must be solved at each time step If these equations are solved by Gaussian elimination then the amount of work required is 0mQN where m is the half bandwidth and N is the total number of equations and unknowns For Poisson s equation on an n by n grid the half bandwidth is m n 1 and the total number of equations is N n 1 for a total of 0n 1 operations For certain problems such 269 as Poisson s equation on a rectangular region with Dirichlet boundary conditions there are fast methods for solving the linear equations that require only 0N log N operations These methods make use of the fast Fourier transform Suppose we have a block TST matrix with TST tridiagonal symmetric Toeplitz blocks ST off 75 A i i S i T i i 1255 T S 3 a 5 7 If A comes from the ve point nite difference approximation 127 for Poisson s equation with Dirichlet boundary conditions on an n by my rectangular grid then with the natural ordering of equations and unknowns S and T are 7 1 by 72 1 matrices given by 2hg 2h 1hg 1hg a hquot 2 12 2 lh n 2h n 2h y where hr 171 and by 171 and A consists of an my 1 by my 1 array of such blocks The eigenvalues and eigenvectors of TST matrices are known It was shown in Theorem 1121 that all TST matrices have the same eigenvectors only their eigenvalues differ Let A be the block TST matrix in 1255 and assume that S and T are ml by ml blocks while the entire matrix is an mg by mg array of such blocks Consider the linear system Au 2 f where the vectors u and f can also be divided into mg blocks each of length ml corresponding to the lines in an ml 1 by Mg 1 grid 11 111 um2T f fl fm2T This linear system has the form Til1 5111 Til1 2 f E 1 Wzg 1256 where uo um21 0 Let Q be the orthogonal matrix whose columns q1qm1 are the eigenvectors of S and T Then S and T can be written as S QAQT T QGQT where A and 8 are diagonal matrices of eigenvalues of S and T respectively Multiplying equation 1256 on the left by QT gives 96239 111 1 MQ39 11 WWII1 Q39fz De ning y E Q39rul and g E Q39I fl this becomes 93 1 1 Aw 9w 1 3 1 quot2 Let the entries in block E of y be denoted 3151 ymhl and similarly for those in each block of g Look at the equations for the jth entry in each block of y 91351 1 My 91371 911 3 1 mz 270 Note that these equations decouple from those for all other entries If for each j 1 ml we de ne Srj 2 3151 yjmmyr and likewise gj 99451 gjsm239r then we have ml independent tridiagonal systems each involving mg unknowns A 694 6 J yjgj g1ml 694 694 A The work to solve these tridiagonal systems is 0mlmg operations This gives us the vectors 3391 Srml which can be rearranged to obtain the vectors yl yma These were obtained from 111 um2 by multiplying by QT Thus the only remaining task is to recover the solution u from y 111 2 9w E 1 ng 1257 This requires mg matrix vector multiplications where the matrices are ml by ml Ordinarily the amount of work required would be 0mgm f operations A similar set of matrix vector multipli cations was required to compute g from f gl 2 Q39l fl E 1 ng It turns out that because of the special form of the matrix Q these matrix vector multiplications can each be performed in time 0ml logg ml using the FFT Thus the total work for performing these multiplications by Q and QT and the major part of the work in solving the original linear system is 0mgml logg ml This is almost the optimal order Omgml 1251 The Fast Fourier Transform To compute the entries of u in 1257 one must compute sums of the form m1 m1 2 wk u x39l 7 sin 7 391 k1 m ch jgilQleM imlJrl 9271 11 01 y y 1 Such sums are called discrete sine transforms and they are closely related to and can be computed from the discrete Fourier transform DFT N l EC ZeW Cij k01N 1 1258 10 Before discussing how to compute the DFT let us brie y describe the Fourier transform in general When one measures a signal f such as the sound made by pushing a button on a touch tone telephone over an interval of time or space one may ask the question what frequencies comprise the signal and what are their magnitudes To answer this question one evaluates the Fourier transform of f Fw e quotftdt co ltwltco 271 A nice feature of the Fourier transform is that the original function f can be recovered from F using a very similar formula for the inverse transform f x27T3C eM Fwdw colttltco The precise de nition of the Fourier transform differs 0 the literature quot the factor 1 m in the de nition of F is not present and then a factor of 127139 is required in the inversion formula Sometimes the de nition of F involves the factor e iwquot and then the inversion formula must have the factor emquot Usually this difference in notation doesn t matter because one computes the Fourier transform performs some operation on the transformed data and then transforms back In practice one cannot sample a signal continuously in time and so one samples at say N points t9 jAt 39 01N 1 and obtains a vector of N function values f1 ftj 39 0 1 N 1 Knowing f at only N points one cannot expect to compute F at all frequencies 11 Instead one might hope to evaluate F at N frequencies wk 2 27rkNAt k 0 1 N 1 The formula for F cuC involves an integral which might be approximated by a sum involving the measured values ff 1 x r 1 quot1 At quot1 N Fwk ewkl t dt z E Z ewk39 39f At E Z ezmk ij oc 10 10 With this in mind the discrete Fourier transform DFT is de ned as in 1258 It is an approximation modulo a factor At 27139 to the continuous Fourier transform at the frequencies wk Like the continuous Fourier transform the DFT can be inverted using a similar formula f1 239 1 N l 7 Z e WWNFk 39 0 1 1 1259 lg0 Evaluating the DFT 1258 or the inverse DFT 1259 would appear to require 0N2 operations One must compute N values F0 FN1 and each one requires a sum of N terms It turns out however that this evaluation can be done with only ON log N operations using the fast Fourier transform FFT algorithm of Cooley and Tukey The ideas of the FFT were implicit in a number of earlier papers such as one by Danielson and Lanczos It was not until the 1960 s however that the complete FFT algorithm was developed and implemented on a computer The key observation in computing the DFT of length N is to note that it can be expressed in terms of two DFT s of length N 2 one involving the even terms and one involving the odd terms in 1258 To this end assume that N is even and de ne in elmN We can write N2 1 N2 1 FlC Z 2W129k7 Vf2j Z 27r1291lc3Vf2j1 10 10 N2 1 N2 1 Z e2m1cN2 fgj wlc Z e2m1cN2 f21 10 10 Fgmwgoz k01N 1 1260 where F a denotes the DFT of the even numbered data points and F 0 the DFT of the odd numbered data points Note also that the DFT is periodic with period equal to the number of data points Thus F and F have period N2 FISH p 0 1 N2 1 Therefore once the N 2 entries in a period of F a and F 0 have been computed we can obtain the length N DFT by combining the two according to 1260 If the N coef cients 217 k 01N 1 have already been computed then the work to do this is 2N operations we must multiply the coef cients wk by the entries Flgo and add to the entries Fly k 0 1 N 1 This process can be repeated For example to compute the length N 2 DFT F a assuming now that N is a multiple of 4 we can write N l l Np Flge Z 62712jlcN2f1j Z 2m211cm392f4j2 10 10 NM l Nl39ll Z e2m1cN ofdj w2c Z eyzm1cN o 1H2 10 10 Fgwumng k01N2 1 where F 8E denotes the DFT of the even even data g and F E0 the DFT of the even odd data fq g Since the length N4 DFT s Fm and Fm are periodic with period N4 once the N 4 entries in a period of these two transforms are known we can combine them according to the above formula to obtain F a This requires 2N2 operations if the powers of w have already been computed and the same number of operations is required to compute F 0 for a total of 2N operations Assuming that N is a power of 2 this process can be repeated until we reach the length 1 DFT s The DFT of length 1 is the identity so it requires no work to compute these There are logQN stages in this process and each stage requires 2N operations to combine the DFT s of length 27 to obtain those of length 2H1 This gives a total amount of work that is approximately 2N log N The only remaining question is how to keep track of which entry in the original vector of input data corresponds to each length 1 transform eg F mi 2 f This looks like a bookkeeping nightmare But actually it is not so dif cult Reverse the sequence of e s and o s assign 0 to e and 1 to o and you will have the binary representation of the index of the original data point Do you see why this works It is because at each stage we separate the data according to the next bit from 273 the right The even data at stage one consists of those elements that have a 0 in the rightmost position of their binary index while the odd data consists of those elements that have a 1 in the rightmost bit of their index The even even data points have a 0 in the rightmost two bits of their index while the even odd data points have a 0 in the rightmost bit and a 1 in the next bit from the right etc To make the bookkeeping really easy one can initially order the data using bit reversal of the binary index For example if there are 8 data points then one would order them as fuuuu fuuuu f1uu1 f 11uu f2u1u f2u1u f3u11 gt f611u f 11uu f1uu1 f5101 f5101 f611u f3u11 f7111 f7111 Then the length 2 transforms are just linear combinations of neighboring entries the length 4 transforms are linear combinations of neighboring length 2 transforms etc 126 Multigrid Methods Exercises 1 Let A denote the Laplacian operator An um aw Show that a problem of the form An f in Q g on 3939 can be solved as follows First nd any function 2 that satis es A2 f in 9 without regard to the boundary conditions Then solve for w Am 2 0 in Q m g 2 on 3939 Then set u 2 217 That is show that 2 w satis es the original boundary value problem 5 9 Suppose you wish to solve Poisson s equation An f on the unit square with u g on the boundary Using a 4 x 3 grid that is 3 interior mesh points in the a direction and 2 interior mesh points in the 3 direction write in matrix form the system of linear equations that you would need to solve in order to obtain a second order accurate approximate solution in On the course web page is a code to solve Poisson s equation on the unit square with Dirichlet boundary conditions Auf in Q 9 on 3939 where Q is the unit square 01 x 01 Download this code and read through it to see if you understand what it is doing Then try running it with some different values of h It uses the same number of subintervals in the a and 3 directions so hr by h Currently it is set up to solve a problem with fa 7 32 32 and g 1 If you do not know an analytic solution to a problem one way to check the code is to solve a problem on a ne grid and pretend that the result is the exact solution then solve on coarser grids and compare your answers to the ne grid solution However you must be sure to compare solution values corresponding to the same grid points Use this idea to check the accuracy of the code Report the approximate L2 norm or in nity norm of the error for several different grid sizes Based on your results what would you say is the order of accuracy of the method 3 Write a code to solve the heat equation in one spatial dimension uLUEm 035531 tgt0 ua0 uga u0t u1t 0 using forward differences in time and centered differences in space lcH 1c M u 1 k k k 1 At 1 27 l EH l u E l Take uga a1 and go out to time t 1 Plot the solution at time t 0 at time t 1 and at various times along the way Experiment with different mesh sizes and different timesteps and report any observations you make about accuracy and stability 5 05 1 Write a code to solve the wave equation 2 LL a an t 2 0 um 0 ai using the difference method 0lt 1 ut0 2955 u0t u1t 0 MI 1c lc l y 2 Jr 61 lt2ultcgt ultlcgt ultlcgt At H2 9 11 1 1 Take a1 3 0 and a 2 2 and go out to time t 2 Plot the solution at each time step you should see the solution behave like a plucked guitar string Experiment with different values of h and At in order to verify numerically the stability condition h At 3 7 la Turn in plots of your solutions for values of h and At that satisfy the stability conditions and also some plots showing what happens if the stability requirements are not met Consider the Lax Friedrichs method 16 44 2 1W1 l 6 6 k j 2 11 1 1 2h At for the differential equation m awn 0 Carry out a stability analysis to determine values of Ath for which the method is stable This is a written exercise Let f be a vector of length 8 with components f0 f1 f7 Write down and clearly explain the sequence of operations you would use to compute the discrete Fourier transform of f using the FFT technique Your explanation should be suf ciently clear that someone unfamiliar with the FFT could follow your instructions and compute the FFT of any given vector of length 8


Buy Material

Are you sure you want to buy this material for

50 Karma

Buy Material

BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.


You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

Why people love StudySoup

Steve Martinelli UC Los Angeles

"There's no way I would have passed my Organic Chemistry class this semester without the notes and study guides I got from StudySoup."

Allison Fischer University of Alabama

"I signed up to be an Elite Notetaker with 2 of my sorority sisters this semester. We just posted our notes weekly and were each making over $600 per month. I LOVE StudySoup!"

Jim McGreen Ohio University

"Knowing I can count on the Elite Notetaker in my class allows me to focus on what the professor is saying instead of just scribbling notes the whole time and falling behind."

Parker Thompson 500 Startups

"It's a great way for students to improve their educational experience and it seemed like a product that everybody wants, so all the people participating are winning."

Become an Elite Notetaker and start selling your notes online!

Refund Policy


All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email


StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here:

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.