4Finite Difference Schemes for the HeatEquationIn the previous chapter we derived a very powerful analytical method forsolving partial differential equations. By using straightforward techniques,we were able to find an explicit formula for the solution of many partialdifferential equations of parabolic type. By studying these analytical solu-tions, we can learn a lot about the qualitative behavior of such models. Thisqualitative insight will also be useful in understanding more complicatedequations.In this section we will turn our attention to numerical methods for solv-ing parabolic equations. Having spent quite some time on deriving elegantanalytical formulas for solving such equations, you may wonder why weneed numerical methods. There are several reasons. We can summarize themain difficulties of the Fourier method as follows:Nonlinear problems. The Fourier method derived above cannot handlenonlinear equations. Both separation of variables and the principle ofsuperposition will in general fail for such equations. This is quite animportant drawback, since many applications give rise to nonlinearproblems. There is a strong tradition for linearizing such equationsin order to apply linear techniques like the Fourier method discussedabove. But this approximation is in some cases very crude. Nonlinearequations can be handled adequately using finite difference methods.The basic ideas are exactly the same as for linear problems, but dif-ficulties may arise in, for example, solving the nonlinear algebraicequations involved.118 4. Finite Difference Schemes for the Heat EquationVariable coefficients. Even linear problems with variable coefficients maybe hard to solve using the Fourier method. In particular this is thecase for discontinuous coefficients. Variable coefficients are fairly easyto handle with finite difference schemes.Integrals. As mentioned above, some of the integrals involved in com-puting the Fourier coefficients can be difficult or even impossible toevaluate analytically. In such cases, numerical integration is needed.Infinite series. In order to actually graph the Fourier solution of a prob-lem, we have to compute the sum of the series. If the series is infinite,we have to rely on an approximation based on a truncation of the se-ries. Furthermore, except for some trivial cases, the partial sum hasto be computed numerically, i.e. on a computer.We conclude that there are a lot of interesting problems that cannot besolved by the Fourier method. And for most problems that can be solved,we are dependent on some sort of numerical procedure in order to actu-ally graph the solution. These observations clearly motivate the study ofnumerical methods in a more general setting. However, you should not bemisled into believing that numerical methods solve every problem. Thereare a lot of dangers in using these methods. The most important difficultieswill be pointed out here.Although there are a lot of different numerical methods available for solv-ing parabolic equations, we focus on finite difference schemes. The reasonfor this is that they are very simple to understand and easy to generalize forquite complicated problems. Furthermore, they are very easy to implementon a computer.In numerical analysis, much effort is spent on the search for efficientschemes in the sense of optimizing accuracy and minimizing CPU timeand memory requirements.1 No claim is made here that finite differenceschemes are optimal in this respect. Often, the most powerful techniquesfor solving partial differential equations use as much analytical informa-tion as possible. The so-called spectral methods are illustrative examplesof this phenomenon. These methods are carefully constructed in order totake advantage of all the analytical insight we have gained through thedevelopment of Fourier techniques. You may consult the notes by Gottlieband Orszag [12] to read more about this topic.The most popular method in applied areas is probably the finite ele-ment method. Although in many respects it is quite similar to the finiteFor scalar computers, i.e., for computers where you only have access to one singleprocessor, it is fairly easy to rank different schemes with respect to these criteria. Justcompute, either a priori or run-time, the number of arithmetic operations needed for thedifferent schemes. In the era of parallel computing this issue is more complicated. Onsuch machines, the quality of a certain method has to be related to how well it exploitsthe access to numerous processors.14.1 An Explicit Scheme 119difference method, the finite element method is preferable when it comes tocomplicated geometries in several space dimensions. You can find a friendlyintroduction to this topic in the book by C. Johnson [15]. A mathematicallymore advanced approach is given by Brenner and Scott [4], and engineersseem to prefer the book by Zienkiewicz [31]. Further references can be foundin these books.4.1 An Explicit SchemeIn this section we will derive a finite difference approximation of the fol-lowing initial-boundary value problem:ut — uxx for x
0,u(0,t) = u(l,t) = 0, u(a;,0) = f{x).(4.1)The way we derive the finite difference scheme for (4.1) is very similar tothe way we derived a scheme for the two-point boundary value problemin Section 2.2. The basic idea is again to replace the derivatives involvedin (4.1) by finite differences. But for this problem we have to approximateboth the space and the time derivatives.Let n > 1 be a given integer, and define the grid spacing in the x-direction by Aa; = l/(w + 1). The grid points in the a;-direction are givenby Xj = jAx for j = 0,1,... ,n+ 1. Similarly, we define tm = mAt forintegers m > 0, where At denotes the time step. Finally, we let vj1 denotean approximation of u(xj,tm).Before we define the scheme, let us recall that we have the followingapproximations2andqilrr AT f\\ Ooifnr i~\\ -A- nl'r -4- AT f}Ax2These approximations motivate the following scheme:_2—At = JAxi for j = 1mBy using the boundary conditions of (4.1), we put<* = 0 and u™+1 = 02See Project 1.1, page 28.120 4. Finite Difference Schemes for the Heat Equationm+1mj-lj+lFIGURE 4.1. The computational molecule of the explicit scheme.for all m > 0. The scheme is initialized byv° = f(xj) for j = 1,... ,n.Let r = Ai/Az2; then the scheme can be rewritten in a more convenientform^+i=TO™1 + (l-27>™+r^+1, j = l,...,n, m>0. (4.2)When the scheme is written in this form, we observe that the values on thetime level im+i are computed using only the values on the previous timelevel tm. Therefore we refer to this scheme as explicit. This is in contrastto implicit schemes where we have to solve a tridiagonal system of linearequations in order to pass from one time level to the next. Such schemeswill be discussed below.In Fig. 4.1, we have sketched the basic structure of this scheme. We oftenrefer to such illustrations as the computational molecule3 of the scheme.Before we start analyzing the properties of the scheme, we present someexamples illustrating how this scheme works.imates one of the particular solutions derived in Section 3.2 above. Thus,we letf(x) = sin(27ra;),and recall from (3.18) that the exact solution is given byu(x,t) =e\"4w2tsm(27rx).3EXAMPLE 4.1 In the first example we look at how well this scheme approx-Sometimes the term \"stencil\" is used for such illustrations.4.1 An Explicit Scheme 121dx= 0.167, dt=0.0125, r=0.450.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9FIGURE 4.2. The finite difference solution (dashed line) and the Fourier-basedsolutions (solid line) of the heat equation. For the finite difference scheme we useda relatively coarse mesh; Ax = 1/6 and At = 1/80.We choose Ax = 1/6 and At = 1/80 and compute the numerical solutionfor 0 < tm < 1/10. The numerical and analytical solution at t = 1/10 isplotted in Fig. 4.2. As usual we have used piecewise linear interpolationbetween the grid points.In Fig. 4.3 we have used a finer mesh, Ax = 1/20 and At = 1/800. Notethat the approximation is much better in this case.EXAMPLE 4.2 In our next example we consider the following initial datafor the problem (4.1):2x 2(1 -x) x < 1/2,x > 1/2.A formal solution of this problem can be derived using Fourier's method.The solution is(4-3)In Fig. 4.4 we have plotted the Fourier solution given by (4.3) as a func-tion of x for t = 0.1. The series is truncated after 200 terms. We also plotthe numerical solution generated by the scheme (4.2) using Ax = 1/50 andAt = 1/5000, hence r = 1/2. In Fig. 4.5, we have plotted the Fourier solu-tion and the numerical solution again, but we have increased the time stepslightly; At = 0.000201. This gives r = 0.5025 and we observe from theplot that something is wrong; the numerical solution oscillates, whereas the122 4. Finite Difference Schemes for the Heat Equationdx= 0,05, dt=0.00125, r=0.5FIGURE 4.3. The finite difference solution (dashed line) and the Fourier-basedsolutions (solid line) of the heat equation. For the finite difference scheme we usedthe mesh parameters Ax = 1/20 and At = 1/800.analytical solution is smooth and very well behaved. This behavior of thescheme is referred to as an instability problem. Much effort will be investedbelow in deriving precise conditions to avoid such problems.4.2 Fourier Analysis of the Numerical SolutionIn our experiments above, we observed two interesting features of the nu-merical solutions. First we noticed that the scheme may generate very goodapproximations of the analytical solutions. And secondly, we saw that bychanging the grid parameters slightly, severe oscillations appeared in thenumerical solution. It is quite clear that such oscillations are unacceptable.The analytical solution is smooth and well behaved, and we look for anapproximation which shares these properties.Our aim in this section is to understand these observations. We want toknow when oscillations can appear and we want to know how to preventsuch behavior. Furthermore, we want to gain some insight into why thescheme generates very good approximations for proper grid sizes. It is notthe purpose of this section to derive any rigorous error analysis for thediscrete solutions. The problem of determining the convergence rate of thescheme will be addressed later. Here we will use our insight in the analyticalsolution to try to understand in a qualitative way why the scheme works.Our analysis of the numerical method will, step by step, follow the pro-cedure outlined in Section 3.1. By using a discrete version of the Fouriermethod, we will derive an explicit formula for the discrete solution. This4.2 Fourier Analysis of the Numerical Solution 1230.4 0.5 0.6 0.7 0.8 0.9FIGURE 4.4. The numerical (dashed line) and Fourier-based solution (solid line)of the heat equation. For the numerical method we have used r = 1/2.formula enables us to compare the analytical and the numerical solutionsterm by term.In the next section, we will present von Neumann's stability analysis.This is a versatile tool which is commonly used to investigate the stabilityof numerical methods for solving time-dependent partial differential equa-tions. The basis of von Neumann's method is the term-by-term analysismentioned above. Both analytical and numerical methods can be decom-posed into linear combinations of particular solutions. Thus, in order tocompare the solutions, it is sufficient to compare the particular solutions.In this section we will do such a comparison thoroughly in order to preparefor the next section. So if you feel that the present section is a bit lengthyand full of details, we promise you that the von Neumann technique wearrive at in the next section is very simple and powerful.4-2.1 Particular SolutionsThe first step in our discrete Fourier analysis is to derive a family of par-ticular solutions of the following problem:v+1At Ax2with the boundary conditionsTfor j = 1,... ,n, m> 0, (4.4)v^ = 0 and <+1 =0, m > 0. (4.5)The initial data will be taken into account later.We are looking for particular solutions of the problem (4.4) and (4.5). Inthe continuous case, we derived the particular solutions by guessing that124 4. Finite Difference Schemes for the Heat Equation* • a0.3« .; •\" ' \" ' ? »i1 »0.25I*Jio4 0.1 ?7''^SNil ?I \"A,\".^0.20.150.12?dx=0.02dt=O000201,r=0.5025 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.05V0.9 a1FIGURE 4.5. The numerical (dashed line) and Fourier-based solution (solid line)of the heat equation. For the numerical method we have used r = 0.5025. Notethat the numerical solution for this value of r oscillates.the space and time dependency could be separated. Thus, we inserted anansatz of the formu(x,t) = X(x)T(t)into the equation. Exactly the same procedure will be applied in the discretecase. Hence we seek particular solutions of the formj1 = XjTm for j = 1,... , n, m > 0.(4.6)Here X is a vector of n components, independent of m, while {Tm}m>o isa sequence of real numbers. By inserting (4.6) into (4.4), we getXjTm-[-i — XjTm ^ Xj^iTm — 2XjTm + Xj+\\TrnAt(Ax)Since we are looking only for nonzero solutions, we assume that XjTm ^ 0,and thus we obtainT-m+i — Tm _ Xj-i — 2Xj + Xj+iAtTm [AjThe left-hand side only depends on m and the right-hand side only dependson j. Consequently, both expressions must be equal to a common constant,say (— n), and we get the following two difference equations:j-i — 2Xj(Az)Atfor j = l,..(4.7)(4.8)-/iTm for m > 0.4.2 Fourier Analysis of the Numerical Solution We also derive from the boundary condition (4.5) thatXo = Xn+1 = 0. 125(4.9)We first consider the equation (4.8). We define4 Tfco = 1, and consider thedifference equationTm+i = (I - Atn)Tm Some iterations of (4.10),Tm+i = (1 - At n)Tm = (1 - A* /x)2Tm_!... ,clearly indicate that the solution isTm = (l-Atn)m for m>0. (4.11)for m>0. (4.10)This fact is easily verified by induction on m.Next we turn our attention to the problem (4.7) with boundary condition(4.9). In fact this is equivalent to the eigenvalue problem (2.44). Hence, fromLemma 2.9 we obtain that the n eigenvalues //i, /i2, • • • , /•*« are given byfik = sin2 (kit Ax/2) for k=l,...,n (4.12)and the corresponding eigenvectors Xk = {Xk,i,Xk,2, • • • >-Xfc,n) S Kn,k = 1,... , n, have components given byXk,j = sin (kitXj) for j = 1,... ,n.Hence, we obtain particular solutions w™j of the formw%j = (1 - Atnk)m sin (knxj). (4.13)So far we have derived a family of particular solutions {wk}%=i with valuesw™j at the grid point (XJ ,tm). Next, we observe that any linear combinationof particular solutionswhere the 7fcS are scalars, is also a solution of (4.4) and (4.5). This obser-vation corresponds to the second step in Section 3.1. Finally, we determinethe coefficients {-ft} by using the initial conditionv° = f{xj) 4for j = 1,... ,n.We are free to choose any nonzero constant; cf. the similar discussion in the contin-uous case on page 91.126 4. Finite Difference Schemes for the Heat EquationSince uik = Xk at t — 0, we want to determine {7^} such thatk.j = f(xj) k=ifor j = 1,... , n.(4.14)Hence, it follows from (2.47) thatxJ)XkJ for fc=l,...,n.(4.15)There is an alternative procedure to derive the representation of thegeneral solution of the finite difference.scheme (4.4)-(4.5). Let A € M.n'nbe the matrix/ 2 -1 A =(Ax)2019 -1 2 0-10 \\01(4.16)V 0and for m > 0 let vm G I\" be the vectorU 0-12/— (Vx ,V2 ,... ,Vn ).The difference scheme (4.4)-(4.5) can then be equivalently written in theformvm+l = (I-AtA)vm, (4.17)where / is the identity matrix. As a simple consequence of this, we obtainvm = (I~ AtA)mv°. (4.18)Recall from Section 2.4 that the eigenvalue problem (2.44) is equivalent tothe eigenvalue problem for the matrix A. Hence, the vectors Xk introducedabove are eigenvectors for the matrix A with corresponding eigenvalue /ifc.Hence, Xk is also an eigenvector for the matrix (/ — AtA)m with eigenvalue(1 - Atfik)m. Hence, if v° = Xk, we derive from (4.18) thatvf = (1 - Atnk)msin(knxj)is a particular solution of (4.4)-(4.5). This solution corresponds exactly to(4.13) above. In the same way as above, the general solution is obtained bytaking linear combinations of these solutions.4.2 Fourier Analysis of the Numerical Solution 1274-2.2 Comparison of the Analytical and Discrete SolutionWe now have explicit formulas of the solutions for both the continuousproblem (4.1) and the discrete problem (4.2). For ease of reference, werepeat the formulas here. The solution5 of the continuous problem is givenbyoou(x, t) = Y fc=icfce~Afc*sin {kirx), (4.19)where \\k = (A;TT)2 andCk = 2 / f{x) sin {kirx) dx. (4.20)JoWe want to compare this analytical solution with the discrete solution givenbynmv — \\ ^ry n _ At tih)msin (kirx ) (4 21)whereandn7^ = 2 Ax \\~] ffaj) sin {k-KXj).for k — 1,... , n. In order to compare the analytical and numerical solutionat a grid point {xj,tm), we define uj1 — u(xj,tm), i.e.oouf = Y, cke~Xktm sin (kirxj). fc=i(4.22)Our aim is to give a rough, but instructive, argument for the fact thatunder appropriate conditions on the mesh parameters Ax and At. To avoidtechnicalities, we consider a fixed grid point (xj,tm) where tm > i fort > 0 independent of the mesh parameters. Furthermore, we assume thatthe initial function / is smooth and satisfies the boundary conditions, i.e.5 ... still formal.128 4. Finite Difference Schemes for the Heat Equation/(0) = /(I) = 0. Finally, we assume that the mesh parameters At and Axare sufficiently small.In order to compare ujl and v™, we note thatU =Tfc=in ckek=\\ Xktmoo ain(kTfXj) + ^ fc=re+cke~eXktmHere, we want to show thatcke~Xktm sin (kirxj) « 0. fe=n+l(4.23)To do this we make the following observations:• Since / is smooth, it is also bounded, and then the Fourier coefficientsCfc given by (4.20) are bounded6 for all k. Hence there is a finiteconstant c such that \\ck\\ < c for all k.• Obviously, we have |sin(fc7TXj)| < 1.By using these observations, we getktm sin t. Since we haveverified (4.23), it follows thatl ckt~Xktm sin (kirxj). fc=i(4.24)Now we want to compare the finite sums (4.24) and (4.21). Motivated bythe derivation of the solutions, we try to compare the two sums termwise.Thus, we keep k fixed, and we want to comparecke~Xktm sin (ktrxj) The Fourier coefficients actually converge towards zero as k tends to infinity. Thisis a consequence of Bessel's inequality, which is discussed in Chapter 8 below. Here,boundedness is sufficient..2 Fourier Analysis of the Numerical Solution and129-yk(l - At fik)m sinikirxj).Since the sine part here is identical, it remains to compare the Fouriercoefficients ck and jk, and the time-dependent terms e~Xktm and (1 -We start by considering the Fourier coefficients, and note that jk is agood approximation of ck becausen2A:r 2_^ f(xj)sm i=i{kirxj)is the trapezoidal-rule7 approximation off12 / f(x) sin (hnx) dx.JoIn fact, by Exercise 2.1 on page 82, we havefor / sufficiently smooth.4-2.3 Stability ConsiderationsFinally, we must compare the time-dependent terms e~*ktm and (l—At/j,k)m.Before we compare the actual values of these expressions, let us briefly con-sider the magnitudes involved. Since Xktm is positive, we havefor all k = 1,... ,n; it is reasonable to require that alsofor all k = 1,... ,n. From this requirement, it follows that Atfik < 2, orequivalentlyfor all k — 1,... ,n. This is always the case if the condition7The trapezoidal method of numerical integration is discussed in Project 2.1 on page82.130 4. Finite Difference Schemes for the Heat Equationis satisfied.Recall now that in Example 4.2 on page 121 we observed that by notobeying the stability condition (4.25), severe oscillations appeared in oursolution. Now we see the reason; if |1 — At/x^l > 1 for some index k,the term (1 - Ai//fc)m blows up as m becomes large. Since such behaviorcannot appear in the analytical solution, we conclude that the condition(4.25) must be satisfied in order to expect reliable numerical results. Thisstability condition will be derived again later using other techniques.4-2-4 The Accuracy of the ApproximationLet us now return to the accuracy considerations initiated above. Of course,we assume from now on that the mesh sizes are chosen such that (4.25)is satisfied. The remaining problem is to discuss how well the term (1 —At/ik)m approximates the term e^Afc'm. In order to study this question, wesimplify the problem a bit by choosing a fixed time tm, say tm = 1, andwe assume that At = (Ax)2/2. As a consequence, we want to compare thetermsak = e~Afcand0k = (I- Atnk)1/At = (1 - 2 sin2 (fc7TV/At/2))1/At.Obviously, ak is very small for large values of k. The first three terms aregiven by«!« 5.172-10\"5, «2« 7.157- 10\"18, a3 « 2.65 • 10~39.The values of 0k depends both on fc and the mesh parameters. By choosinga relatively fine mesh, say Ax = 1/100, we get/?!« 5.1-10~5, 02 « 6.973 • 10\"18, 03 « 2.333 • MT39.These computations clearly indicate that both ak and 0k become very smallwhen k is large. Furthermore, we observe that 0k seems to approximate akadequately. We will consider this problem a bit more closely. Since both akand 0k are very small for large values of k, it is sufficient to compare themfor small k.4.2 Fourier Analysis of the Numerical Solution 131In order to compare u^ and 0k for small values of k, we start by recallingthatsin (y) =y + O(y3).Thus, we get2 sin2 (A;7r-v/Ai/2) « {kw)2At.Furthermore, we have in general thatfor e sufficiently small. By using these facts we derive0k = (1 - 2 sin2 (fc7TV/Ai/2))1/A*w (1 - (fc7r)2Ai)lM*This shows that also the time dependent term e~Xktm is well approximatedby its discrete counterpart, (1 — Atfik)m-4-2.5 Summary of the ComparisonIn order to summarize our analysis, we go back to the exact solution givenbyooXkt~Xkt sin (kirx),u(x, t) = ^2 cte~fc=iand the representation of the discrete solution given byWe have seen thatn(i) u(xj,tm) « yj Cfce~Afetm sin (kirXj)by truncating the Fourier series. In (i), the Fourier coefficients satisfy(ii) c^ — il -1 nf(x) sin (/CTTO;) dx ss 2Ax V^ /(^j) sin (fcTTXj) = 7^•/o -1132 4. Finite Difference Schemes for the Heat Equationby the trapezoidal rule of numerical integration. Furthermore, if the meshparameters satisfy(Ax)2 ~then(in) e~Xktm « (]by the properties of the eigenvalues.The observations (i), (ii), and (in) imply thatu(Xj,tm) « ^7fc(l \" Atfik)msm(knXj) = vThis explains why we get good approximations for appropriate choices ofAa; and At.We have derived a stability condition (4.25) which has to be satisfiedin order to get well-behaved numerical approximations. Secondly, we haveutilized a discrete version of the Fourier method to show that each sig-nificant term in the analytical solution is well approximated by a similarterm in the discrete solution. Although this analysis is rough and does notsupply a precise error estimate, it explains very well what is going on inour computations. And furthermore, it is useful in order to prepare for thevon Neumann stability analysis that we will present in the next section.The basic idea of this technique is exactly the same as in the present sec-tion; stability of the entire approximation is studied by analyzing particularsolutions. It is important to note that this way of analyzing a scheme ispurely linear; no similar method is found for general nonlinear problems.Therefore, we will come back to other ways of investigating the stability ofnumerical methods later.4.3 Von Neumann's Stability AnalysisIn the section above, we saw that a discrete version of the Fourier methodenabled us to understand important features of a numerical scheme. Thebasic observation is that questions regarding stability and convergence canbe analyzed by comparing the analytical and numerical solutions term byterm. Thus, particular solutions of both problems become very important.In this section we will continue along the same lines. Our aim is to generalizethis way of investigating the stability of a scheme to cover a wider class ofequations and boundary conditions.4.3 Von Neumann's Stability Analysis 1334-3.1 Particular Solutions: Continuous and DiscreteLet us start by recalling some particular solutions of the heat equation,ut = uxx for x e (0,1), t> 0.In the presence of Dirichlet boundary conditions,u(0, i) = u(l,t) — 0, t > 0,the particular solutions are given byFQ = {Tk{t) sin (krrx)}™-!,whereFor Neumann data,ux{0,t)=ux(l,t)=0, the particular solutions are given byt>0,see Section 3.6. And finally, for periodic boundary conditionsu(-l,t)=u(l,t) and ux(-l,t) = ux(l,t), t>0,where the space variable x € (—1,1), we haveFP = FDUFN;cf. Exercise 3.15 on page 111.In order to be able to handle all these particular solutions in a uniformmanner, it is convenient to write them in a slightly different form. We knowfrom calculus thatelx = cos (x) + i sin (x), where i is the imaginary unit. Hence, we havecos (a;) = -(eix +e~lx)andSm(4.26)^~2Te ~6 ''134 4. Finite Difference Schemes for the Heat EquationUsing these formulas, all the functions in the families Fp, F^: and Fp canbe expressed as linear combinations of the following functions8- {JfcWe }fc=_oo-In a similar way, we can argue that the corresponding discrete problemshave a family of particular solutions of the formF 1f/« \\mAkiTXi 1 oowhere ak represents the time dependency of the discrete solutions. In theexplicit scheme for the heat equation, it is given by ak — 1 — At/ik', see(4.21). This term is often referred to as the amplification factor of thescheme.4-3.2 ExamplesThe basic idea of von Neumann's method is to compare the growth ofthe analytical and discrete particular solutions. More precisely, we wantto derive conditions on the mesh parameters Ax and Ai such that thegrowth of the discrete solutions are bounded by the growth of the analyticalsolutions. Let us look at two examples to clarify this procedure.EXAMPLE 4.3 We consider the heat equationut=uxx. By inserting a particular solution of the form(4.27)into (4.27), we getTk(t) = -(kn)2Tk(t).Hence,Tfc(t) = e~{k^2\\ (4.28)where we as usual have defined Tk(0) = 1. The solution of the partialdifferential equation is approximated by the schemeV+1T - VT _v?~i- At ~ 2vT + VJ+IAl? 'Why do we suddenly need complex functions? So far, everything has been real. Donot worry too much about this; the introduction of complex variables here is merely atool to simplify our calculations. We can do this directly by using the sine and cosinefunctions, or we can handle both at once by using the complex exponentials.84.3 Von Neumann's Stability Analysis By inserting a particular solution of the form135we get(ak)m+1 At - (ak)m 6lk _ eik**>-i - ~Since Xj = jAx, this impliesAiCOS yf£7Tl\\X }- 1(Ax)4 n. 22(Ax) °Hencewe haveak = 1 -2m4At .2(Ax)2CTTAX/2)SinceTk,given by (4.28), satisfies\\Tk(t)\\ < 1for allk,we also require thatfor all k. As above, cf. (4.25), this inequality holds if the following conditionis satisfied:Thus, for both Dirichlet, Neumann, and periodic boundary conditions, thiscondition has to be satisfied in order to get reasonable numerical results.•EXAMPLE 4.4 Let us apply the procedure to the following equation:ut = uxx + u. Again, by inserting the particular solution(4.30)136 4. Finite Difference Schemes for the Heat Equationinto the equation (4.30), we getrk(t) = (1 - HenceTk(t) = and thus, for all k, we have\\Tk(t)\\ 0. (4.31)e^-W*(kn)2)Tk(t).The solution of equation (4.30) is approximated by the schemeAt By insertingAx2 J 'we getfDiknA.xAt or(Ax)2(4.32)Based on the bound (4.31) for the analytical solution, it is reasonable torequire that the numerical solution satisfiesfor all k. Suppose that the usual conditionAt(Ax)2is satisfied. Then by (4.32) we get0, 137(4.34)which holds for all y > — 1. The proof of this fact is left to the reader inExercise 4.27. The conclusion of this example is that, again, the condition(4.33) must be satisfied in order to get stable results. •We summarize our discussion so far by stating that a numerical solutionis said to be stable in the sense of von Neumann if the growth of the dis-crete particular solutions can be bounded by the growth of the continuousparticular solutions. More precisely, if we letT(t)=m&x\\Tk(t)\\,kthen we say that the scheme is stable in the sense of von Neumann ifmax|(afe)m| 0.4-3.3 A Nonlinear ProblemAs mentioned above, the method of von Neumann is only valid for lin-ear problems with constant coefficients. That, of course, is a major draw-back, because linear problems with constant coefficients are about the onlyproblems we can solve analytically. So, under the circumstances where webadly need numerical methods, e.g. for nonlinear problems and problemswith variable coefficients, the method of von Neumann cannot be applieddirectly. However, in practical computations, the method of von Neumannis applied far beyond the family of problems where we have actually shownthat the method works. It is frequently applied to both linear problems withvariable coefficients and to nonlinear problems. The procedure is roughlyas follows: Given a nonlinear problem, we linearize the equation and freeze9the coefficients by considering the problem locally. In this manner, we derivea linear problem with constant coefficients. For this problem, the methodof von Neumann can be applied, and a stability condition can be derived.Certainly, this condition will depend on the frozen coefficients involved.The trick is then to choose a conservative time step, covering all possiblevalues of the frozen coefficient. More precisely, we try to find coefficientvalues that lead to the most restrictive time step. In some cases, this canbe quite difficult, since we do not know any bounds for the coefficients. One Freezing the coefficient means to approximate the coefficients by a constant. Ofcourse, this can only be valid locally, and thus freezing of coefficients often leads to afamily of problems with different constant coefficients.9138 4. Finite Difference Schemes for the Heat Equationpractical solution of this problem, is to introduce a variable time step, andupdate the bounds on the coefficients at each time step.10Let us illustrate this procedure by an example.EXAMPLE 4.5 Consider the following nonlinear heat equation,ut = (a(u)ux)x u(0,t) = ti(l,t)=O,for x 6 (0,1), t > 0,u{x,0) = f{x),whereAn explicit finite difference scheme for this problem is given by,Ax2At ^J~^, (4.35)where am+1,2 = [a{vJ+ i) + a(vJl))/2- The derivation of this scheme will beconsidered in Exercise 4.20.Consider this problem locally, i.e. close to some fixed location (xo,to).If u is smooth, we can approximate the function a(u) by a constant valuea0 = a(u(xa, to)) close to (x0, t0). This approximation leads to the equationut = aouxx,and the associated scheme— — = «0 ~2 f°r J^l)---)^! m > 0.The particular solutions of the linearized equation are given byTk(t)elk™,whereTk{t) = e\"^2^1satisfies the usual bound|Tfc(t)| < 1, 10t > 0,The technique of linearizing the equation and freezing the coefficients can be carriedout in order to analyze amazingly complicated problems. An excellent introduction tothis procedure is given in the book by Kreiss and Lorenz [17]. A thorough discussion ofthe practical use of von Neumann's method for complicated problems can be found inthe book by Godunov and Ryabenkii [10].4.3 Von Neumann's Stability Analysis 139for all k since ao > 1. Consequently, we require that the particular solutionsof the corresponding finite difference scheme,satisfy the bounda-kfor all k. By inserting the discrete particular solution into the scheme, wefind that,4aoAt 2( , .Thus we require that the mesh parameters satisfy the bounda0At/(Ax)2 < 1/2.This heuristic argument indicates that for mesh parameters satisfying thisbound, the scheme is stable, at least locally. In order to derive a globalcondition, we observe thata{u) 1+3u2 =< 3for all u. Thus any frozen coefficient ao is less than 3, and consequentlythe most restrictive requirement on the time step is given byAt<{-^. (4.36)In Fig. 4.6, we have plotted the numerical solution when the initial datais given by f(x) = sin(37ra), using the mesh parameters At — 0.00005and Ax = 0.02. We observe that the solution seems to be well behaved.The reader is encouraged to do further experiments using this scheme; seeExercise 4.19.We will return to the problem considered in this example later11 andactually prove that the requirement (4.36) is a sufficient condition for thenumerical solutions to be stable. As indicated earlier, we will derive tech-niques that are more suitable for nonlinear problems. However, the presentexample indicates that by doing some rough arguments, the von Neumannmethod can be used to derive reasonable time-step requirements even fornonlinear problems. Of course, time steps derived by this procedure mustbe applied with great care. •Further examples of stability analysis based on von Neumann's methodwill be given in the exercises.nSee Section 6.3.2 on page 190.140 4. Finite Difference Schemes for the Heat EquationFIGURE 4.6. The finite difference solution of the nonlinear heat equation usingAx = 0.02 and At = 0.00005.4.4 An Implicit SchemeWe have studied one particular numerical method for solving the heat equa-tion. The method is given by (4.2) and is referred to as explicit. Obviously,the scheme is very simple to implement on a computer, and we have seenthat it has some nice properties. Further properties will be derived below.However, the explicit method suffers from one major drawback; it requiresvery small time steps due to the stability condition.Let us look a bit closer on the consequences of the stability conditions(4.25), i.e.At< 1/2.(4.37)Suppose we want to compute a numerical solution of the heat equation attime t = 1, and that accuracy requirements force us to choose a fairly finemesh, say n = 100. Then by the stability condition we must haveA£<120402'Since we want the solution at time tu = 1, we must take M = 20 402time steps. Refining the mesh by choosing n = 1000, we have to computeM = 2 004 002 time steps. Clearly, even this very simple problem can puteven our most powerful modern computers under strain. Of course, byturning our interest towards two or three space dimensions, this situationbecomes even worse.This unfortunate feature of the explicit scheme motivates us to searchfor alternatives with higher computational efficiency. Indeed there are a4.4 An Implicit Scheme 141lot of methods available. In this section, we will present the simplest andprobably most popular method: the standard implicit scheme. We do notwant to go too deep into the discussion of different methods here; the topicis far too large and the interested reader may consult e.g. Thomee [27] andreferences given there for further discussions.Before we start presenting the implicit scheme, let us remind ourselvesof the basic difference between explicit and implicit schemes. We statedabove, on page 120, that a scheme is called explicit if the solution at onetime step can be computed directly from the solution at the previous timestep. On the other hand, we call the scheme implicit if the solution on thenext time level is obtained by solving a system of equations.We want to derive an implicit scheme for the following equation:ut — uxx for x g (0,1), t > 0,u(0,t)=0, u(l,i)=0, u(x,0) = f(x).Borrowing the notation from the explicit scheme, we apply the followingapproximations:i + _L Art ~ u(x,t + At)-u(x,t)ut(x,t + At) «Atand. . . u(x — Ax, t + At) — 2u(x, t + At) + u(x + Ax, t + At)Kuxx(x, t + At) « -i ! '- -^—2 >- ^ '- '-.AxThis leads to the following scheme:(4.38)Vi+vT=Vi+ito~ ~C2v+Vj++1 for j=i>--->'n m^°-The computational molecule of this scheme is depicted in Fig. 4.7.The boundary conditions of (4.38) imply thatv'S1 = 0 and u™+1 = 0for all m > 0, and the initial condition givesv°j= f{xj) for j = 1,... ,n.In order to write this scheme in a more convenient form, we introduce thevector vm g W1 with components vm = (n™,... , i>™)T. Then we observethat the scheme can be written as(J + AtA)vm+l =vm, m> 0. (4.39)142 4. Finite Difference Schemes for the Heat Equationm+1mj-lj+lFIGURE 4.7. The computational molecule of the implicit scheme.where / 6 K™-™ is the identity matrix, and where the matrix A g M.n'nisgiven by (4.16) above, i.e./ 2 -1 A =-1 2 '•• 0 ... -1 \"-. ••• ••. 0 \\:o(AXyo (4.40): V 0 '•• ... -1 2 0 -1 -12 /In order to compute numerical solutions based on this scheme, we haveto solve linear systems of the form (4.39) at each time step. Hence, it isimportant to verify that the matrix (/ + AtA) is nonsingular such thatmvm+i js uniqUeiy determined by v. In order to prove this, we use theproperties of the matrix A derived in Lemma 2.9 on page 70.Lemma 4.1 The matrix (I + AtA) is symmetric and positive definite forall mesh parameters.Proof: The matrix (I + AtA) is obviously symmetric, since A is symmetric.Furthermore, the eigenvalues of (I + AtA) are of the form 1 + Atfj,, where /icorresponds to eigenvalues of A. However, the eigenvalues of A, which aregiven by (4.12), are all positive. Therefore, all the eigenvalues of (I + AtA)are positive, and hence this matrix is positive definite. •Since (I + AtA) is symmetric and positive definite, it follows from Propo-sition 2.4 that the system (4.39) has a unique solution that can be computedusing the Gaussian elimination procedure given in Algorithm 2.1 on page53. From a computational point of view, it is important to note that thecoefficient matrix (I + AtA) in (4.39) does not change in time. This obser-4.4 An Implicit Scheme 143dx=0.02, dt=0.000201, r=0.5025FIGURE 4.8. The numerical solution (dashed line) computed by the implicitscheme and Fourier-based solution (solid line) of the heat equation. For the nu-merical method we have used r = 0.5025.vation can be used to reduce the total amount of computational effort inthe scheme; see Exercise 2.13 on page 75.Let us see how the implicit scheme works.EXAMPLE 4.6 In Example 4.2 we observed that for one choice of grid pa-rameters, the explicit scheme produces very good approximations. However,by increasing the timestep slightly, severe oscillations appear. Now we wantto see how the implicit scheme handles this situation.In Fig. 4.8 we have plotted the analytical solution, computed as in Ex-ample 4.2, and the numerical solution provided by the implicit scheme.The grid parameters are given by Ax = 0.02 and At = 0.000201, thusr = At/(Ax)2 = 0.5025. We recall that these parameters gave an oscilla-toric solution using the explicit scheme. From the figure, we observe thatthe numerical solution computed by the implicit scheme is very well be-haved.This observation leads us to believe that reliable solutions can be com-puted by this scheme without obeying the stability condition (4.37). Let usgo one step further and choose At = Ax = 0.02, which gives r = 50. Theresults are given in Fig. 4.9, and we see that the numerical solution is stillnice and smooth.4-4-1 Stability AnalysisWe observed in the previous example that the stability condition derived forthe explicit scheme seems unnecessary for getting reasonable results usingthe implicit scheme. Let us look closer at this phenomenon by applying144 4. Finite Difference Schemes for the Heat Equation0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1FIGURE 4.9. The numerical (dashed line) and Fourier-based solution (solid line)of the heat equation. For the numerical method we have used r = 50.the von Neumann method derived above. Recall first that the particularsolutions of the heat equation (4.38) are given bywhereBy inserting a particular solution of the forminto the implicit scheme (4.39), we getr, 1 £> — iknAx o _i_ ^ihnAxkAt(Ax)2 (Ax)Consequently,- sin2 (knAx/2).11+(£JI sin2 (fcTrAx/2)'By observing that\\Tk(t)\\ <4.5 Numerical Stability by Energy Arguments for all k, we require|(afc)ro| < 1- Since145(4-41)1 +where all the /x^s are strictly positive, it follows that the requirement (4.41)is fulfilled for all mesh parameters. This explain why oscillations did notappear in the computations above.Numerical methods that are well behaved for any choice of grid param-eters are referred to as unconditionally stable. We have just seen that theimplicit scheme deserves this label. The explicit scheme is analogously re-ferred to as conditionally stable.It should be noted here that although we are able to compute approx-imations using arbitrarily long time steps, the issue of accuracy has notyet been discussed. Obviously, by choosing At very large, the accuracy ofthe computation is poor. In numerical analysis this issue is a topic of livelydiscussion; should implicit or explicit schemes be used? The problem is ofcourse how to compute good approximations using as little CPU time andmemory resources as possible.124.5 Numerical Stability by Energy ArgumentsIn Section 3.7 we introduced energy arguments in order to derive a stabilityproperty for the solution of the heat equation. A similar analysis can alsofrequently be performed for finite difference solutions. It is possible to derivecertain properties of the solution of the finite difference method withoutknowing the solution in detail. Here we shall illustrate these techniques bystudying the solution of the explicit finite difference method (4.4) applied tothe initial and boundary value problem (4.1). Recall that if r = At/(Ax)2,this difference scheme has the formv?+1=v? + r(v?_1-2v? + v?+1), j = l,...,n, m > 0, (4.42)with boundary conditions< = 0. (4.43) You might think that with the extremely powerful computers of today, such consid-erations are less important than earlier. This is not the case. We always strive for higheraccuracy and more complicated models far beyond the capacity of any known computer.146 4. Finite Difference Schemes for the Heat EquationFurthermore, we will assume throughout this section that the stability con-dition (4.25) holds, i.e.1 - 2r > 0. (4.44)Along the same lines as in Section 3.7, we introduce, for each time levelm > 0, the discrete energyE = Ax^ivJ1)2, mn(4.45)and we are interested in the dynamics of this scalar variable. More preciselywe want to show that E decreases with time, i.e.Em+1 QInstead of computing the time derivative of the energy, as we did in Section3.7, we consider the corresponding time difference,Em+1 _£m=Al3 == AxJT{v™+x+vf){v™+l - v?)= rAx J2(v?+l + \"DKli - 2< + v?+1) (4.47)3 = 13=1 3=1where we have used the difference scheme (4.42).We consider each of the three parts on the right-hand side separately.Observe first that from the boundary condition (4.43) and by summationby parts (cf. (2.31) on page 60) we obtain?>£.! - 2vf3 = 1 3 = 1Furthermore, by applying the difference scheme (4.42) and summation byparts once more,n3=1 )=1- z 2-i^ 3 > + 3 = 1 /r 2-J^V3 + l 3 = 1V3 I •4.5 Numerical Stability by Energy Arguments Finally, we use the inequality13ab< -(a2 + b2)- 2to obtain147E«r+i«-i+«r+i) ^ X>r+1)2+l^uf+«+i)2))3=1 3=13 = 1Collecting these three inequalities, it follows from (4.47) thatEm+iwhere we have used the stability assumption (4.44). Hence,+1 -Em) <0,and by (4.44) this implies the desired inequality (4.46).We summarize the result of the discussion above:Theorem 4.1 Let {vj1} be a solution of the finite difference scheme (4-42)-(4-43) and let the corresponding energy {Em} be given by (4-45). If thestability condition (4-25) holds, then {Em} is nonincreasing with respect tom.Hence, we have seen that the stability condition (4.25), or (4.44), impliesthat the explicit difference scheme admits an estimate which is similar tothe estimate (3.60) for the continuous problem. As for the continuous prob-lem, this can be used to estimate the difference between two solutions, withdifferent initial data. This follows since the difference of two solutions ofthe finite difference scheme is a new finite difference solution. We thereforeobtainCorollary 4.1 Assume that the stability condition (4-25) holds and let{vj1} and {wj1} be two solutions of the finite difference scheme (4-42)-(443). Then, for all m > 0,£>f - wff < Ax3=1 3=113See Exercise 4.24.148 4. Finite Difference Schemes for the Heat EquationThe interpretation of this result is that the difference scheme is a stabledynamical system in the sense that an error in the initial data boundsthe error in the corresponding solutions. Energy arguments can also beperformed for the implicit scheme (4.39). This is discussed in Exercise 4.25below.4.6 ExercisesEXERCISE 4.1 Verify, by direct calculation, that the discrete functions{wn} given by (4.13) are solutions of (4.4)-(4.5).EXERCISE 4.2 Implement the scheme (4.2) for the heat equation and inves-tigate the performance of the method by comparing the numerical resultswith the analytical solution given by(a) Example 3.1 on page 92.(b) Example 3.2 on page 93.(c) Example 3.4 on page 97.pare the numerical solutions provided by the explicit and the implicitschemes.EXERCISE 4.4 In this exercise we want to study the rate of convergence ofEXERCISE 4.3 Repeat Exercise 4.2 using the implicit scheme (4.39). Com-the explicit scheme by doing numerical experiments. Define the error byeA(tm)= j=0,... ,n+lmax \\U(XJ, tm) - v™\\.JWe want to estimate the rate of convergence for the scheme at time t = 1/10for the problem considered in Exercise 4.2 (a).(a) Estimate, using numerical experiments, a such that eA(l/10) = O((At)a)for At = (Ax)2/2.(b) Repeat the experiments in (a) using At = (Ax)2/6.(c) Try to explain the difference in the rate of convergence encounteredin the two cases above. Hint: Consider the truncation error discussedin Exercise 4.15 below.4.6 Exercises EXERCISE 4.5 Consider the following initial-boundary value problem149ut = uxx for x e (0,1), t > 0,u(0,i) = «<(*), u(l,t)=ur(t),u(x,0) = f(x).Here ue(t), ur(t), and /(a;) are bounded functions satisfying U((0) = /(0)and ur(0) = /(I).(a) Derive an explicit scheme for this problem.(b) Derive an implicit scheme for this problem and show that the linearsystem that arises can be solved by Gaussian elimination.EXERCISE 4.6 Derive an explicit scheme for the following Neumann prob-lem:ut — uxx for x e (0,1), t > 0,ux{0,t) = ux(l,t) = 0,u(x,0) = f(x).Use the analytical solution given in Example 3.5 on page 101 to check thequality of your approximations.of the problem. Compare the numerical solutions provided by the explicitand the implicit schemes.EXERCISE 4.8 Consider the problemEXERCISE 4.7 Repeat Exercise 4.6 by deriving an implicit approximationut=auxx u(x,0) = f{x).for xe (-£,£), t> 0,where a, b and I, a > 0 are given constants.(a) Derive an explicit scheme.(b) Derive an implicit scheme.(c) Find the exact solution when a = 2, £ = IT, a = —TT, b = TT, andf(x) = x + sin(3x).(d) Implement the schemes derived in (a) and (6) and compare the resultswith the analytical solution derived in (c).150 4. Finite Difference Schemes for the Heat EquationEXERCISE 4.9 Consider the problemut = 4uxx-10u + q(x,t) for u(£ut) = a(t), u(£2,t) = b(t),u(x,0) = f{x),x e {£i,£2), t > 0,where £2 > £\\ are given constants, and a(t),b(t), and q(x, t) are givenfunctions.(a) Derive an explicit scheme.(b) Derive an implicit scheme.(c) Suppose*i = -2, ^2 = 3, a{t) = et-2, andq(x,t) = lie* + 10a;.Show thatu(x,t) = e* +a;is an exact solution of the problem.(d) Implement the schemes derived in (o) and (b) and compare the resultswith the analytical solution derived in (c).b(t)=et + 3, f(x) = l+x,EXERCISE 4.10 Consider the problemut = (a{x,t)ux)x+c{x,t)ux+q(x,t) u(eltt) = a(t), u(£2,t)=b(t),u(x,0)=f(x),for xe(£i,£2), t > 0,where £2 > £\\ are given constants, and a(t),b(t),a(x,t),c(x,t), are given smooth functions.(a) Derive an explicit scheme.(b) Derive an implicit scheme.and q(x,t)4.6 Exercises EXERCISE 4.11 Consider the problem151ut = (a(x,t)ux)x+c(x,t)ux+q(x,t)u ux{eut) =a(t), ux(e2,t) = b(t),u(x,0) = f(x),for x E {li,l2), t > 0,where £2 > ^1 are given constants, and a(t),b(t),a(x,t),c(x,t), are given functions.(a) Derive an explicit scheme.(b) Derive an implicit scheme.and q(x,t)EXERCISE 4.12 Consider the equationUt = auxx,with Dirichlet boundary conditions. Here a > 0 is a given constant. Wedefine an explicit schemem+l V_• J m — VA J mVA Ar 1 —^Aa;for j = l,...,n, m > 0,and an implicit schemevm+l _ yrn yrn+l(a) Derive a stability condition for the explicit scheme using the vonNeumann method.(b) Show that the implicit scheme is unconditionally stable in the senseof von Neumann.EXERCISE 4.13 Consider the equationut = uxx,with Neumann-type boundary conditions and the following explicit schemev^-j > = vri-2^+vT+1 for j=hn: m>0(a) Use the Taylor series to explain the derivation of this scheme.152 4. Finite Difference Schemes for the Heat Equation(b) For what mesh sizes is this scheme stable in the sense of von Neu-mann?EXERCISE 4.14 Consider the equationu = utxx — 9WJwith Dirichlet-type boundary conditions. Derive an explicit and an implicitscheme for this equation. Use the von Neumann method to investigate thestability of the methods.EXERCISE 4.15 In Section 2.3.5 we introduced the concept of truncationerror for a finite difference approximation for a two-point boundary valueproblem. Here we shall discuss a similar concept for difference approxima-tions of the heat equation.Observe that the scheme (4.4) can be written in the formj-(vm+l -vm) + Avm = 0,where A G M\"'n is given by (4.16).The truncation vector rm G K\" is given byrm = l-(um+1-um) where um G W1 is given by u\"1 = u(xj,tm) problem (4.1).+ Aum,for a solution of the continuous(a) Show that under suitable smoothness assumptions on the solution u,|rf | = O(At) + O((Ax)2). (4.48)In rest of this exercise we study a more general difference scheme of theformj-(vm+1 - vm) + 6{Avm+l) + (1 - 6)Avm = 0, (4.49)where 6 G [0,1] is a parameter. Note that if 6 = 0, this corresponds to theexplicit scheme (4.4), while if 6 = 1, it corresponds to the implicit schemestudied in Chapter 4.4 above.(b) Sketch the computational molecule for the scheme when 6 G (0,1).(c) Show that for all 6 G [0,1] the estimate (4.48) holds, and that thechoice 6=1/2 leads to an improved estimate of the form|rf | = 0({At)2)+0({Ax)2).(Hint: Consider Taylor expansions at the point (XJ, (tm+i + tm)/2).)4.6 Exercises 153EXERCISE 4.16 Motivated by the result of Exercise 4.15, we study, in thisexercise, the difference scheme 4.49 with 0 = 1/2. This difference schemeis usually referred to as the Crank-Nicholson scheme. In component formthe scheme is given byvrn+l _ vrn ^At \"~ 2 ^ (AX)2 + {Axf )for j = 1, 2,... , n and m > 0.(a) Show that this implicit scheme is unconditionally stable in the senseof von Neumann.(b) Discuss how the vectors vm+1 € Wl can be computed from vm.(c) Show that the solution of the Crank-Nicholson scheme for the initial-boundary value problem (4.1) admits the representationVJ =fc=lwhere O(JU) = (l — ^M) (l + ^M) and7fe = 2Ax'=1 v(d) Show that the amplification factor of the difference scheme, a(n),satisfiesK^-e-^'l =O((Ai)3).How does this result relate to the corresponding result for the explicitscheme (4.4)? Compare your result with the conclusions you derivedin Exercise 4.15.(e) Implement the Crank-Nicholson scheme. Choose the initial functionf(x) as in Example 4.2 and try to verify that the scheme is uncondi-tionally stable by varying the parameter r = ,^2.EXERCISE 4.17 Consider the general scheme (4.49). Use the von Neumannmethod to discuss the stability for any 9 e [0,1].EXERCISE 4.18 Consider the equationut +cux = uxx,with Dirichlet-type boundary conditions. Here c > 0 is given constant.(a) Show that this problem has a family of particular solutions of theform— (ikTtc+(kir)2)t ikvx1 4. Finite Difference Schemes for the Heat Equation(b) Show thatTk(t) = e-C^+W)'satisfies the bound14\\Tk(t)\\0,for all k.Derive stability conditions for the following numerical methods by applyingthe von Neumann method.(c)m+l mJ V± ,At c2Aa; Ax2A* (e)At (f)vrn+l _ vm +CAx Ax2+ C-Ax Ax2ymAt (g)vrn+l _ „„ ++ CC A^ Ax2^rn+1 _ 2vm+l =+ ^rn+1ym+l __ ^rn+1 ~^ +C A^ 'A?tion (see (4.36)) based on some rough considerations. The purpose of thisexercise is to perform a numerical study of the quality of this condition.Consider the initial condition f(x) = 100a;(l — x)\\x — 1/2|, and run thescheme (4.35) with several different grids. For this initial function, is thecondition (4.36) sufficient in order to guarantee well-behaved numericalsolutions?14EXERCISE 4.19 In Example 4.5 on page 138, we derived a stability condi- Recall that for a complex number z = x + iy, the absolute value, or the modulus, isgiven by \\z\\ = (x2 + y2)l/2. It is also useful to note that \\ei$\\ = 1 for any 6 eR.4.6 Exercises scheme for the following problem:ut = (a(u)ux)x for x € (0,1), t > 0,u(0,t) = u(l,t) = 0,where a{u) is a given strictly positive and smooth function,(a) Put v = a{u)ux and justify the following approximations:ut(x,t)vx(x,t)(b) Show thatv(x + Ax/2,t) « \\{a(u{x + Ax,t))\\ + a(u{x,t))))u{x(c) Use these approximations to derive the scheme15,,,771+1 v155EXERCISE 4.20 The purpose of this exercise is to derive a finite differenceu(x,t + At) - u(x,t)Atv(x + Ax/2, t) - v(x - Ax/2, t)Axj ~j vTO _ anm (vp1 Vj + l/2\\j + l ~j Vitm\\ )~ nm aAt ~ Ax2j-l/2VjV(ni\"1where a™+l/2 = {a{v™+l) + a{v?))/2.EXERCISE 4.21 Consider the initial-boundary value problem in Exercise4.20. Derive an implicit scheme and investigate the stability of the methodby the technique discussed in Example 4.5 on page 138.EXERCISE 4.22 Consider the nonlinear heat equationih = (ae(u)ux)x for x £ (0,1), t > 0,u(0,t) =u(l,t) = 0,u(x,0) = sin(37rx).Here ae(u) = 2 + ecos (u), where e is a small parameter, |e| m Vj + l/2(j + l — i\\m\\ Vj ) ~j-\\/2\\j ar\\m V (i}m _ ?>m V^At \"\" a{vf))/2.Ax2j-1)where af+1/2 = (a(v™+1) + (b) Leta(u) — u, a(t) = t, b(t) = 1 + t, and f(x) — x.Show that u(x, t) = x + t is an exact solution of this problem.(c) Show, by induction, that the explicit scheme gives the exact solutionat each grid point, i.e. show that vj1 = Xj + tm, for any grid sizes.(d) Compute the numerical solution at t = 1 using the scheme imple-mented in (a) for the problem defined in (6). Try the following gridparameters:- n = 4 and At = 1/65.- n = 30 and At = 1/10.Discuss your observations in light of the result in (c).(e) From the numerical results obtained in (d), it is clear that some kindof stability condition is needed. Use the procedure discussed in Exam-ple 4.5 on page 138 to derive a stability condition for this problem.Run some numerical experiments with mesh parameters satisfyingthis condition. Are the numerical solutions well-behaved if the con-dition on the mesh parameters is satisfied?4.6 Exercises 157EXERCISE 4.24 Use the fact that (a - b)2 > 0 to show thatab< ha2 + b2)for all real numbers a and b.given by (4.40) is positive definite to show that the implicit differencescheme (4.39) satisfies an estimate of the form (4.46) for all mesh parame-ters.EXERCISE 4.26 Similar to the discussion in Chapter 2 we introduce a dis-EXERCISE 4.25 Use an energy argument and the fact that the matrix Acrete inner product which is an analog of the inner product (•, •) for con-™ we define16tinuous functions. For a vectors v, w € K.Hence, this inner product is just the ordinary Euclidean inner productmultiplied by the scaling factor Ax. We recall from Chapter 2 that thisinner product arises naturally when the vector v has the interpretation of adiscrete function defined on the grid points Xj = jAx. We let || • ||A denotethe corresponding norm, i.e.As above we let Xk = (XkA,Xk,2,... vectors with components given by, Xk,n) 6l\ k = 1,... , n, be theXk,j = sin (knXj) for j = 1,... , n.Recall that these vectors are orthogonal with respect to the inner product(•, -)A and that ||Xfc||A = 1/2 (see Lemma 2.30).(a) Explain why any vector v 6 W1 can be written in the formnfc=lwhereck 16=2{v,Xk}A.In Chapter 2 we used h to indicate the spacing in the x-variable, and hence we usedthis subscript to indicate the corresponding discrete inner product. Here, where we havetwo grid parameters Aa; and At, we use A for the same purpose.158 4. Finite Difference Schemes for the Heat Equation(b) Show that(c) Let {i>m}m>o be a sequence of vectors generated by the finite differ-ence scheme (4.39). As above, in Exercise 4.25, let= \\\\vm\\\\l.Show thatwhere fi\\ = ^5 sin2(7rAa;/2).(d) Explain whylim Ern = 0,m—^00and compare this result with what you derived in Exercise 4.25 above.EXERCISE 4.27 Show that(1 + y)m < emyfor all m > 0 and y > — 1.
