FINITE DIFFERENCE METHODS FOR SOLVING DIFFERENTIAL …

1 FINITE DIFFERENCE METHODSFORSOLVING DIFFERENTIAL EQUATIONSI-Liang ChernDepartment of MathematicsNational Taiwan UniversityMay 16, 20132 Contents1 FINITE DIFFERENCE Approximation .. Basic Numerical METHODS for Ordinary DIFFERENTIAL equations .. Runge-Kutta METHODS .. Multistep METHODS .. Linear DIFFERENCE equation .. Stability analysis .. Zero Stability ..182 FINITE DIFFERENCE METHODS for Linear Parabolic FINITE DIFFERENCE METHODS for the Heat Equation .. Some discretization METHODS .. Stability and Convergence for the Forward Euler method .. von Neumann Analysis .. Energy method .. Stability Analysis for Montone Operators Entropy Estimates .. Entropy estimate for backward Euler method .. Existence Theory .. Existence via forward Euler method .. A Sharper Energy Estimate for backward Euler method.

2 Relaxation of errors .. Boundary Conditions .. Dirichlet boundary condition .. Neumann boundary condition .. The discrete Laplacian and its inversion .. Dirichlet boundary condition .. 383 FINITE DIFFERENCE METHODS for Linear elliptic Discrete Laplacian in two dimensions .. Discretization METHODS .. The 9-point discrete Laplacian .. Stability of the discrete Laplacian .. Fourier method .. Energy method .. 444 FINITE DIFFERENCE Theory For Linear Hyperbolic A review of smooth theory of linear hyperbolic equations .. Linear advection equation .. Linear systems of hyperbolic equations .. FINITE DIFFERENCE METHODS for linear advection equation .. Design techniques .. Courant-Friedrichs-Levy condition .. Consistency and Truncation Errors .. Lax s equivalence theorem.

3 Stability analysis .. Modified equation .. FINITE DIFFERENCE schemes for linear hyperbolic systemwith constant coefficients .. Some design techniques .. Stability analysis .. FINITE DIFFERENCE METHODS for linear systems with variable coefficients .. 645 Scalar Conservation Physical models .. Traffic flow model .. Burgers equation .. Two phase flow .. Basic theory .. Riemann problem .. Entropy conditions .. Rieman problem for nonconvex fluxes .. Uniqueness and Existence .. 746 FINITE DIFFERENCE Schemes For Scalar Conservation Major problems .. Conservative schemes .. Entropy and Monotone schemes .. 807 FINITE DIFFERENCE METHODS for Hyperbolic Conservation Flux splitting METHODS .. Total Variation Diminishing (TVD) .. Other Examples for ( ).

4 Extensions .. High Order Godunov METHODS .. Piecewise-constant reconstruction .. piecewise-linear reconstruction .. Multidimension .. Splitting Method .. Unsplitting METHODS .. 998 Systems of Hyperbolic Conservation General Theory .. Rarefaction Waves .. Shock Waves .. Contact Discontinuity (Linear Wave) .. Physical Examples .. Gas dynamics .. Riemann Problem of Gas Dynamics .. 1109 Kinetic Theory and Kinetic Kinetic Theory of Gases .. Kinetic scheme .. 1172 CONTENTSC hapter 1 IntroductionThe goal of this course is to provide numerical analysis background for FINITE DIFFERENCE methodsfor SOLVING partial DIFFERENTIAL equations . The focuses are the stability and convergence theory. Thepartial DIFFERENTIAL equations to be discussed include parabolic equations , elliptic equations , hyperbolic conservation FINITE DIFFERENCE ApproximationOur goal is to appriximate DIFFERENTIAL operators by FINITE DIFFERENCE operators.

5 How to performapproximation? What is the error so produced? We shall assume the underlying functionu:R Ris smooth. Let us define the following FINITE DIFFERENCE operators: Forward DIFFERENCE :D+u(x) :=u(x+h) u(x)h, Backward DIFFERENCE :D u(x) :=u(x) u(x h)h, Centered DIFFERENCE :D0u(x) :=u(x+h) u(x h) ,his called the mesh size. By Taylor expansion, we can get u (x) =D+u(x) +O(h), u (x) =D u(x) +O(h), u (x) =D0u(x) +O(h2).We can also approximateu (x)by other FINITE DIFFERENCE operators with higher order errors. Forexample,u (x) =D3u(x) +O(h3),34 CHAPTER 1. INTRODUCTION whereD3u(x) =16h(2u(x+h) + 3u(x) 6u(x h) +u(x 2h)).These formulae can be derived by performing Taylor expansion ofuatx. For instance, we expandu(x+h) =u(x) +u (x)h+h22u (x) +h33!u (x) + u(x h) =u(x) u (x)h+h22u (x) h33!u (x) + .Substracting these two equations yieldsu(x+h) u(x h) = 2u (x)h+2h33!

6 U (x) + .This givesu (x) =D0u(x) h23!u (x) + =D0u(x) +O(h2).In general, we can derive FINITE DIFFERENCE approximation foru(k)at specific pointxby the valuesofuat some nearby stencil pointsxj,j= 0,..,nwithn k. That is,u(k)(x) =nXj=0cju(xj) +O(hp k+1)for somepas larger as possible. Here, the mesh sizehdenotesmax{|xi xj|}. As we shall seelater that we can choosep=n. To find the coefficientscj,j= 0,..,n, we make Taylor expansionofu(xj)about the pointx:u(xj) =pXi=01i!(xj x)iu(i)(x) +O(hp+1).We plug this expansion formula ofu(xj)into the FINITE DIFFERENCE approximation formula foru(k)(x):u(k)(x) =nXj=0cjpXi=01i!(xj x)iu(i)(x) +O(hp k+1).Comparing both sides, we obtainnXj=01i!(xj x)ihicj= 1ifi=k0otherwise ,fori= 0,..,pThere arep+ 1equations here, it is natural to choosep=nto match then+ 1unknowns. This is an nVandermonde system. It is nonsingular ifxiare different.

7 The matlab code fdcoeffV(k,xbar,x)can be used to compute these coefficients. Reference: Randy LeVeque s book and his Matlab BASIC NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS5In the case of uniform grid, using central FINITE differencing, we can get high order approxima-tion. For instance, letxj=jh, wherej Z. Letuj=u(xj). Thenu (0) =u1 u 1h+O(h2)u (0) =u1 2u0+u 1h2+O(h2)u(3)=12h3(u2 2u1+ 2u0 2u 1+u 2) +O(h2). Considerxi=ih,i= 0,..,n. Let x=xm. Find the coefficientsciforu(k)( x)and thecoefficient of the leading truncation error for the following cases: k= 1,n= 2,3,m= 0,1,2,3. k= 2,n= 2,m= 0,1, Basic Numerical METHODS for Ordinary DIFFERENTIAL EquationsThe basic assumption to design numerical algorithm for ordinary DIFFERENTIAL equations is the smooth-ness of the solutions, which is in general valid provided thecoefficients are also smooth.

8 Basicdesignning techniques include numerical interpolation, numerical integration, and FINITE methodEuler method is the simplest numerical integrator for ODEy =f(t,y)( )is discretized byyn+1=yn+kf(tn,yn).( )Here,kis time step size of the discretization. This method is called the forward Euler method. Itsimply replacedy/dt(tn)by the forward FINITE DIFFERENCE (yn+1 yn)/k. By Taylor expansion, thelocal truncation error is n:=y (tn) y(tn+1) y(tn)k=O(k)Leten:=yn y(tn)be the true C1and suppose the solutiony =f(t,y)withy(0) =y0exists on[0,T]. Then the Euler method converges at anyt [0,T]. In fact, the true errorenhas thefollowing estimate:|en| e t O(k) 0,asn .( )Here, = max| f/ y|andnk= 1. the regularity of the solution, we havey C2[0,T]andy(tn+1) =y(tn) +kf(tn,y(tn)) +k n.( )Taking DIFFERENCE of ( ) and ( ), we obtain|en+1| |en|+k|f(tn,yn) f(tn,y(tn))|+k| n| (1 +k )|en|+k| n|.

9 Where|f(t,x) f(t,y)| |x y|.The FINITE DIFFERENCE inequality has a fundamental solutionGn= (1 + k)n, which is positiveprovidedkis small. Multiplying above equation by(1 + k) n 1, we obtain|em+1|G m 1 |em|G m+kG m 1| m|.Summing inmfromm= 0ton 1, we get|en| n 1Xm=0Gn m 1k| m| n 1Xm=0 GmO(k2)=Gn 1G 1O(k2) Gn O(k) e t O(k),wheret=nkand we have used(1 + k)n e The theorem states that the numerical method converges in[0,T]as long as the solutions ofthe ODE The error isO(k)if the solution is inC2[0,T].3. The proof above relies on the existence and smoothness of the solution. However, one canalso use this approach to prove the local existence theorem by showning the approximatesolutions generated by the Euler method form a Cauchy Euler methodIn many applications, the system is relaxed to a stable solution in a very short period of time. Forinstance, considery = y y.

10 The corresponding solutiony(t) yast O( ). In the above forward Euler method, practically,we should require1 +k BASIC NUMERICAL METHODS FOR ORDINARY DIFFERENTIAL EQUATIONS7in order to haveGnremain bounded. Here, is the Lipschitz constant. In the present case, = 1/ .If is very small, the the above forward Euler method will require very smallkand lead to inefficientcomputation. In general, forward Euler method is inefficient (require smallk) ifmax f(t,y) y >> the case f/ y >>1, we have no choice to resolve details. We have to take a very , if f/ y <0, say for example,y = ywith >>1. then the backward Euler methodis recommended:yn+1=yn+kf(tn+1,yn+1).The error satisfies|en+1| en k|en+1|+O(k2)The corresponding fundamental solution isGn:= (1 + k) n. Notice that the error satisfies|en| n 1Xm=0(1 + k) mO(k2) (1 + k) n+1 kO(k2) e T O(k).There is no restriction on the size frog methodWe integratey =f(t,y)fromtn 1totn+1:y(tn+1) y(tn 1) =Ztn+1tn 1f( ,y( ))d.