Numerical Solution Of Partial Di Erential Equations
Numerical solution of partial differentialequationsDr. Louise Olsen-KettleThe University of QueenslandSchool of Earth SciencesCentre for Geoscience ComputingE–mail: l.kettle1@uq.edu.auWeb: nKettleISBN: 978-1-74272-149-1
AcknowledgementsSpecial thanks to Cinnamon Eliot who helped typeset these lecture notes inLATEX.Note to readerThis document and code for the examples can be downloaded .Please note if any of the links to code are not working please contact myemail address and I can send you the code.
Contents1Overview of PDEs1.1 Classification of PDEs . . . . .1.1.1 Elliptic . . . . . . . . . .1.1.2 Hyperbolic . . . . . . .1.1.3 Parabolic . . . . . . . .1.2 Implicit Vs Explicit Methods to1.3 Well-posed and ill-posed PDEs. . . . . . . . . . . . . . . . . . . . . . . . .Solve PDEs. . . . . . .INumerical solution of parabolic equations2Explicit methods for 1-D heat or diffusion equation2.1 Analytic solution: Separation of variables . . . . . . . . .2.2 Numerical solution of 1-D heat equation . . . . . . . . . .2.2.1 Difference Approximations for Derivative Terms inPDEs . . . . . . . . . . . . . . . . . . . . . . . . .2.2.2 Numerical solution of 1-D heat equation using thefinite difference method . . . . . . . . . . . . . . .2.2.3 Explicit Forward Euler method . . . . . . . . . . .2.2.4 Stability criteria for forward Euler method . . . . .2.3 Method of lines . . . . . . . . . . . . . . . . . . . . . . .2.3.1 Example . . . . . . . . . . . . . . . . . . . . . . . .3Implicit methods for 1-D heat equation3.1 Implicit Backward Euler Method for 1-D heat equation . .3.1.1 Numerical implementation of the Implicit BackwardEuler Method . . . . . . . . . . . . . . . . . . . . .3.1.2 Dirichlet boundary conditions . . . . . . . . . . . .3.1.3 Mixed boundary conditions . . . . . . . . . . . . .3.2 Crank-Nicolson Scheme . . . . . . . . . . . . . . . . . . .2.99991010101213. . 13. . 15. . 15.161720212123. . 23.24242426
45IIIterative methods4.1 Jacobi method . . . . . . . . . . . . . . . .4.1.1 Example . . . . . . . . . . . . . . . .4.1.2 Applying the Jacobi method . . . . .4.2 Gauss-Seidel Method . . . . . . . . . . . . .4.2.1 Example: using Gauss-Seidel methodmatrix equation . . . . . . . . . . . .4.3 Relaxation Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .to solve a. . . . . . . . . . . . . . 31. . . . 322-D Finite Difference5.1 2-D Poisson’s equation . . . . . . . . . . . . . . . . . . . . .5.2 2-D Heat (or Diffusion) Equation . . . . . . . . . . . . . . .5.2.1 Alternating Direct/Implicit method for the 2-D heatequation . . . . . . . . . . . . . . . . . . . . . . . .5.3 Cylindrical and spherical polar co-ordinates . . . . . . . . .5.3.1 Example: Temperature around a nuclear waste rod .34. 34. 38. 39. 40. 42Numerical solution of hyperbolic equations6Analytical solutions to the 1-D6.1 1-D Wave equation . . . . . .6.2 d’Alembert’s solution . . . . .6.3 Separation of variables . . . .7Flux conservative problems7.1 Flux Conservative Equation . . . . . . . . . . . . . . . .7.2 Stability analysis of numerical solutions of the first orderflux conservative or 1-D advection equation . . . . . . .7.3 Forward Time Centred Space (FTCS) . . . . . . . . . .7.3.1 von Neumann stability analysis of FTCS method7.4 Lax Method . . . . . . . . . . . . . . . . . . . . . . . . .7.4.1 von Neumann Stability Analysis of Lax Method7.5 Courant Condition . . . . . . . . . . . . . . . . . . . . .7.6 von Neumann Stability Analysis For Wave Equation . .7.6.1 Lax method . . . . . . . . . . . . . . . . . . . . .7.7 Other sources of error . . . . . . . . . . . . . . . . . . .7.7.1 Phase Errors (through dispersion) . . . . . . . . .7.7.2 Dispersion in the numerical solution of the 1-Dadvection equation using the Lax method . . . .7.7.3 Error due to nonlinear terms . . . . . . . . . . . .3Wave. . . . . . . . . .282929303147equation48. . . . . . . . . . . . . . 48. . . . . . . . . . . . . . 48. . . . . . . . . . . . . . 4951. . . 51.51525253535455555656. . . 57. . . 59
7.7.4 Aliasing error . . . . . . . . . . . . . . . . . . . . . . . 5989Numerical Solution of 1-D and 2-D Wave Equation8.1 Explicit Central Difference for 1-D Wave Equation . . . . .8.1.1 Example: plucking a string . . . . . . . . . . . . . .8.1.2 1-D Wave Equation with Friction . . . . . . . . . . .8.2 2-D Wave Equation . . . . . . . . . . . . . . . . . . . . . . .8.2.1 Example: vibrations of a thin elastic membrane fixedat its walls . . . . . . . . . . . . . . . . . . . . . . . .8.2.2 Examples of wave equation . . . . . . . . . . . . . .Finite element method9.1 An introduction to the Finite Element Method . . . . . .9.2 Comparing FEM solution to FD solution for our example9.2.1 FD solution . . . . . . . . . . . . . . . . . . . . . .9.3 2-D Finite Element Method . . . . . . . . . . . . . . . . .9.3.1 2-D “hat functions” . . . . . . . . . . . . . . . . . .9.3.2 Example: 2-D Finite Element Method using eScriptfor elastic wave propagation from a point source. .6161616568. 68. 71.737377777879. . 8010 Spectral methods8310.1 An introduction to spectral methods . . . . . . . . . . . . . . 8310.1.1 Example 1: Comparing the accuracy of solutions of avariable speed wave equation with either the spectralor finite difference method . . . . . . . . . . . . . . . . 8410.1.2 Example 2 Comparing spectral and finite differencemethods with constant wave speed conditions andinitial conditions of a non-smooth pulse . . . . . . . . . 87IIINonlinear partial differential equations11 Shock wave11.1 Analytical solution: Method of characteristics . . . . . . . . .11.1.1 Example 1: Using method of characteristics to solvethe linear 1-D advection equation . . . . . . . . . . . .11.1.2 Example 2: Using method of characteristics to solvethe nonlinear inviscid Burger’s equation . . . . . . . .11.2 Numerical Solution for nonlinear Burger’s Equation . . . . . .11.2.1 Example I: Finite difference solution with Lax Method11.2.2 Example II: Solution using Method of Lines . . . . . .48990909191949597
11.2.3 Example III: Solution using Spectral Method . . . . . . 9712 Korteweg-de Vries Equation12.1 Solitons . . . . . . . . . . . . . . . . . . . . . .12.2 Analytical solution . . . . . . . . . . . . . . . .12.3 Numerical solution of KdV Equation . . . . . .12.3.1 Solving directly with Spectral Method .12.3.2 Modifying Uxxx term causing instabilitiesspectral method . . . . . . . . . . . . . .12.3.3 Interacting Solitons . . . . . . . . . . . .5. . . . . . . . . . . . . . . . .in direct. . . . . . . . .100. 100. 100. 103. 104. . . 105. . . 106
List of Figures2.1 Initial conditions in (a) and matlab solution using ForwardEuler method for temperature distribution along rod withtime in (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.1 Initial conditions in (a) and matlab solution using BackwardEuler method for temperature distribution along rod withtime in (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.1 Plot of residual using the Gauss-Seidel method after eachiteration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325.1 3-D Cylindrical Co-ordinates . . . . . . . . . . . . . . . . . . . 415.2 3-D Spherical Polar Co-ordinates . . . . . . . . . . . . . . . . 415.3 Initial conditions in (a) and matlab solution using BackwardEuler method for temperature distribution near nuclear rodat different time intervals in (b) . . . . . . . . . . . . . . . . . 467.1 Solution at t 1 using the Lax method with different timesteps, (a) t x/2c where dispersion is present but thepulse matches the analytical solution for the speed of thewave, (b) t x/c where no dispersion is present andnumerical solution matches analytical solution exactly, and(c) t 1.001 x/c where the Courant condition is not metand solution is becoming unstable. . . . . . . . . . . . . . . . 597.2 Aliasing error occurs when the mesh spacing x is too largeto represent the smallest wavelength λ1 and misinterprets itas a longer wavelength oscillation λ2 . . . . . . . . . . . . . . 608.1 Initial conditions in (a) and matlab solution using explicitcentral difference method for 1D wave equation in (b) . . . . 648.2 D’Alembert’s solution in (a) and error using numericalmatlab solution using explicit central difference method for1D wave equation in (b) . . . . . . . . . . . . . . . . . . . . . 656
8.3 Initial conditions in (a) and matlab solution using explicitcentral difference method for 1D wave equation with frictionin (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 689.1 FEM mesh with triangles . . . . . . . . . . . . . . .9.2 FEM mesh with triangles . . . . . . . . . . . . . . .9.3 2D hat function(φj (xj , yj ) 1, φj (xi , yl ) 0 if i 6 j and j 6 l) . .9.4 Plot of Euclidean normal of the displacement at t 0point source using eScript. . . . . . . . . . . . . . . . . . . 74. . . . . 78. . . . . 80for a. . . . . 8110.1 Numerical solution for 1D advection equation with initialconditions of a smooth Gaussian pulse with variable wavespeed using the spectral method in (a) and finite differencemethod in (b) . . . . . . . . . . . . . . . . . . . . . . . . . . 8710.2 Numerical solution for 1D advection equation with initialconditions of a box pulse with a constant wave speed usingthe spectral method in (a) and finite difference method in (b) 8811.1 The analytical solution U (x, t) f (x U t) is plotted toshow how shock and rarefaction develop for this example .11.2 Initial conditions in (a) and solution for nonlinear Buger’sequation using the Lax method in (b) . . . . . . . . . . . .11.3 Initial conditions in (a) and solution for nonlinear Buger’sequation using the method of lines in (b) . . . . . . . . . .11.4 Initial conditions in (a) and solution for nonlinear Buger’sequation using the spectral method in (b) . . . . . . . . .q. . 95. . 96. . 98. . 99Ax). . . . .12.1 The initial conditions U (x, 0) f (x) Asech2 ( 1212.2 Initial conditions in (a) and solution for nonlinear KdVequation using the direct spectral method in (b) . . . . . . .12.3 Initial conditions and final solution after one period in (a)and solution for nonlinear KdV equation using a modifiedspectral method in (b) . . . . . . . . . . . . . . . . . . . . .12.4 The initial conditions U (x, 0) f (x) is plotted to show the 2solitions and their speeds . . . . . . . . . . . . . . . . . . . .12.5 Initial conditions and final solution after one period in (a)and solution for nonlinear KdV equation for two interactingsolitons in (b) . . . . . . . . . . . . . . . . . . . . . . . . . .7. 102. 105. 106. 107. 107
Numerical solution of parabolicand hyperbolic PDEsThis document and code for the examples can be downloaded .Please contact me (l.kettle1@uq.edu.au) if you are unable to downloadthe code and I can send it to you.References: Applied Numerical Methods for Engineers using Matlab and C, R. J. Schillingand S. L. Harris. Computational Physics Problem solving with computers, R.H. Landauand M. L. Páez. An Introduction to Computational Physics, T. Pang. Numerical Recipes in Fortran (2nd Ed.), W. H. Press et al. Introduction to Partial Differential Equations with Matlab, J. M. Cooper. Numerical solution of partial differential equations, K. W. Morton andD. F. Mayers. Spectral methods in Matlab, L. N. Trefethen8
Chapter 1Overview of PDEs1.1Classification of PDEsThe classification of PDEs is important for the numerical solution you choose.A(x, y)Uxx 2B(x, y)Uxy C(x, y)Uyy F (x, y, Ux , Uy , U )1.1.1EllipticAC B 2For example, Laplace’s equation:Uxx Uyy 0A C 1, B 01.1.2HyperbolicAC B 2For example the 1-D wave equation:Uxx A 1, C 1/c2 , B 091Uttc2
1.1.3ParabolicAC B 2For example, the heat or diffusion EquationUt βUxxA 1, B C 01.2Implicit Vs Explicit Methods to Solve PDEsExplicit Methods: possible to solve (at a point) directly for all unknown values in thefinite difference scheme. stable only for certain time step sizes (or possibly never stable!). Stability can be checked using Fourier or von Neumann analysis. Timestep size governed by Courant condition for wave equation.Implicit Methods: there is no explicit formula at each point, only a set of simultaneousequations which must be solved over the whole grid. Implicit methods are stable for all step sizes.1.3Well-posed and ill-posed PDEsThe heat equation is well-posed Ut Uxx . However the backwards heatequation is ill-posed : Ut Uxx at high frequencies this blows up!In order to demonstrate this we let U (x, t) an (t) sin(nx)then:Uxx an (t)n2 sin(nx),U U t {z xx}andUt ȧn (t) sin(nx) ȧn (t) sin(nx) an (t)n2 sin(nx)Heat Equation10
2tȧn an n2 an (t) an (0)e nFor the heat equation the transient part of the solution decays and this hasstable numerical solutions.U t Uxx{z ȧn (t) sin(nx) an (t)n2 sin(nx)}Backwards Heat Equation2tȧn an n2 an (t) an (0)enFor the backwards heat equation the transient part of the solution blows upand the numerical solution would fail! In general it is difficult or impossibleto obtain numerical solutions for ill-posed PDEs.11
Part INumerical solution of parabolicequations12
Chapter 2Explicit methods for 1-D heator diffusion equationWe will focus on the heat or diffusion equation for the next few chapters.This is an example of a parabolic equation.2.1Analytic solution: Separation of variablesFirst we will derive an analtical solution to the 1-D heat equation. Considerthe temperature U (x, t) in a bar where the temperature is governed by theheat equation, Ut βUxx . The ends of the bar are cooled to 0 C and theinitial temperature of the bar is 100 C.U (0, t) 0 C- U (L, t) 0 C6U (x, 0) 100 CWe want to solve Ut βUxx using separation of variables. We assumethat the solution can be written as the product of a function of x and afunction of t, ie. U (x, t) X(x)T (t) then: 2X TX βT βUxx divide by XT t x2X 00 (x)2(2.1) λ {z}X(x) {z }constantonlyfunction of x onlyUt 1 T 0 (t)β T (t) {z }function of t13
The only way the LHS and RHS of equation 2.1 can be a function of t andx respectively is if they are both equal to a constant which we define to be λ2 for convenience.2T 0 λ2 βT 0 T e λ βtX 00 λ2 X 0 X A sin λx B cos λx Use boundary conditions U (t, 0) 0U (t, L) 0t X(0) 0 B X(L) 0 A sin λnπλ λn ,n 1, 2, . . .LU (x, t) X(x)T (t) XAn sin(n 1nπx λ2 βt)eL2e λ βt is a transient solution and decays in time to boundary conditions.Use initial conditions U (0, x) 100 C to find An :U (0, x) T0 XAn sin (nπx/L)n 1Use orthogonality: 0L sin( nπx) sin( mπx)dx δnmLLand cos(mπ) 1 0, for m 0, 2, 4, . . .and cos(mπ) 1 2, for m 1, 3, 5, . . .R Am T0 [ L/mπ(cos(mπ) 1)] 2LT0 mπform 1,3,5,. . .14
2.2Numerical solution of 1-D heat equation2.2.1Difference Approximations for Derivative Termsin PDEsWe consider U (x, t) for 0 x a,0 t TDiscretise time and spatial variable x: t tk k t,T,m x a,n 10 k m xj j x,0 j n 1Let Ujk U (xj , tk )Consider Taylor series expansion for Ujk 1 :Ujk 1 Ujk t Ujk t2 2 Ujk O( t3 ) t2 t2(2.2)If we only consider O( t) terms in equation 2.2 then we arrive at the forwarddifference in time approximation for Ut : UjkUjk 1 Ujk O( t) t tWe can also derive a higher order approximation for Ut if we consider theTaylor series expansion for Ujk 1 as well:Ujk 1 Ujk t Ujk t2 2 Ujk O( t3 ) t2 t2(2.3) UjkUjk 1 Ujk 12.2 2.3 O( t2 ) leap-frog (or centred difference) in time. t2 tThis gives higher order accuracy than forward difference.15
We can also perform similar manipulations to arrive at approximationsfor the second derivative Utt :2.2 2.3 2Ujk 2 UjkUjk 1 Ujk 1 2Ujk O( t2 ) central difference t2 t2The finite difference method makes use of the above approximations tosolve PDEs numerically.2.2.2Numerical solution of 1-D heat equation usingthe finite difference methodUt βUxxInitial conditionsU (0, x) f (x)Types of boundary conditionsNeumann boundary conditionsUx (t, 0) g1 (t)Ux (t, a) g2 (t)Dirichlet boundary conditionsU (t, 0) g1 (t)U (t, a) g2 (t)Mixed boundary conditionsU (t, 0) g1 (t)Ux (t, a) g2 (t)16
2.2.3Explicit Forward Euler methodor FTCS (Foward Time Centred Space)We want to solve the 1-D heat equation:Ut βUxx .(2.4)We solve this PDE for points on a grid using the finite difference methodwhere we discretise in x and t for 0 x a and 0 t T :xn 1 axn.x3x2x1x0.t0 t1 t2 t3tm 1 tm TWe discretise in time with time step: t T /m and in space with gridspacing: x a/(n 1), and let tk k t where 0 k m and xj j xwhere 0 j n 1.Let Ujk U (xj , tk ) then the finite difference approximations for equation(2.4) are given by:Ujk 1 Ujk U (xj , tk ) O( t), Forward Euler method for time derivative, t tkkUj 1 2Ujk Uj 1 2 U (xj , tk ) O( x2 ), x2 x2Central difference method for spatial derivative.Our discretised PDE (equation (2.4)) becomes: Ujk 1 Ujkβ kkk U 2U Ujj 1 ,j 12 t x kk Uj 1 (1 2s)Ujk ,or Ujk 1 s Uj 1where s β t. x217
Ujk 1 is the solution for the temperature at the next time step.Suppose we have initial conditions U (x, 0) Uj0 f (xj ), and mixedboundary conditions:Dirichlet boundary conditions at x 0: U (0, t) U0k g1 (tk ),k Un 1Neumann boundary conditions at x a: Ux (a, t) g2 (tk ), xNumerical implementation of Explicit Forward Euler methodSolving equation (2.4): Ut βUxx with:initial conditions: Uj0 f (xj ) U (x, t 0),Dirichlet boundary conditions at x 0: U (x 0, t) U0k g1 (tk ) andkNeumann boundary conditions at x a: U (a, t)/ x Un 1/ x g2 (tk ).To solve using the Neumann boundary condition we need an extra step:kk Un 1 UnkUn 1 g2 (tk ), x xkor Un 1 xg2 (tk ) Unk .(2.5)We can write out the matrix system of equations we will solve numericallyfor the temperature U . Suppose we use 5 grid points x0 , x1 , x2 , x3 , x4 xn 1 ,ie. n 3 in this example:x1x0 0x2x3x4 xn 1 aWe let:U1k k k at time tk . U2 , solution for temperature vector UkU3 kU kThe boundary conditions give U0k U (x 0, tk ) and Un 1 U4k U (x a, tk ). kkWe can rewrite Ujk 1 s Uj 1 Uj 1 (1 2s)Ujk in matrix form:U1k 11 2ss0U1ksU0k k 1 k s1 2ss U2 U2 0 (2.6)0s1 2sU3ksU4kU3k 1 k 1U 18
Using boundary conditions: U0k g1 (tk ) and U4k xg2 (tk ) U3k equation(2.6) becomes: k 1 U U1k 1U2k 1U3k 1 1 2sss1 2s0s {zA Dirichlet bc sg1 (tk )0s xg2 (tk )z } { {z}0s1 {z s} Neumann bc }U1kU2k U3k Neumannbc } {z b k 1or U AU b k(2.7)The term (1 s) in the matrix A above and the term (s xg2 (tk )) invector b above are from the Neumann boundary condition given using theapproximation in equation (2.5).Matlab code for Explicit Forward Euler methodThe matlab code for this example is ForwardEuler.m.Solving equation (2.4): Ut βUxx with 0 x 1, β 10 5 , and0 t 12, 000, using m 600 time steps, and n 39 for 41 grid points inx.Initial conditions: U (x, t 0) 2x sin(2πx) 1,and Dirichlet boundary conditions at x 0: U (x 0, t) 1 and Neumannboundary conditions at x 1: U (x 1, t)/ x 2.With these boundary conditions we can check that the numerical solutionapproximates the steady state solution U (x, t) 2x 1 as t . Yoursolution using the finite difference code can also be checked using Matlab’sPDE solver (using pdex1).Figure 2.1 shows the initial conditions in (a) and matlab solution for temperature distribution along rod with time in (b). The numeri
Numerical Recipes in Fortran (2nd Ed.), W. H. Press et al. Introduction to Partial Di erential Equations with Matlab, J. M. Cooper. Numerical solution of partial di erential equations, K. W. Morton and D. F. Mayers. Spectral methods in Matlab, L. N. Trefethen 8
The main objective of the thesis is to develop the numerical solution of partial difierential equations, partial integro-difierential equations with a weakly singular kernel, time-fractional partial difierential equations and time-fractional integro partial difierential equations. The numerical solutions of these PDEs have been obtained .
the numerical solution of second-order optimization methods. Next step development of Numerical Multilinear Algebra for the statistical analysis of multi-way data, the numerical solution of partial di erential equations arising from tensor elds, the numerical solution of higher-order optimization methods.
Numerical solution of partial differential equations in science and engineering. "A Wiley-Interscience publication." Includes index. 1. Science—Mathematics. 2. Engineering. mathematics. 3. Differential equations, Partial— Numerical solutions. I. Pinder, George Francis, 1942- II. Title. Q172.L36 515.3'53 81-16491 ISBN 0-471-09866-3 AACR2
No analytical solution exists or it is difficult to obtain it. Solution of Nonlinear Equations Solution of Linear Equations Curve Fitting Least Squares Interpolation Numerical Integration Numerical Differentiation Solution of Ordinary Differential Equations Solution of Partial Differential Equations
A Numerical Solution to the Partial Di erential Equation of the Risk Based Capital for Guaranteed Minimum Withdrawal Bene t Zhan Zhang, YuFei Hou, Yang Xu June 8, 2015 Abstract This research project is dedicated to implement an algorithm for the partial di erential equation solution to a risk management problem of the GMWB variable annuity .
Indo-German Winter Academy, 2009 3 Need for Numerical Methods for PDE’s Most of the PDEs are non-linear Most of them do not have analytical solutions Difficult to find analytical solution in most cases due to its complexity Even if the analytical solution can be found, computing it takes more time than that needed for numerical solution
Preface to the First Edition The book is designed for use in a graduate program in Numerical Analysis that is structured so as to include a basic introductory course and subsequent more specialized courses. The latter are envisaged to cover such topics as numerical linear algebra, the numerical solution of ordinary and partial differential .
American Chiropractic Board of Radiology Heather Miley, MS, DC, DACBR Examination Coordinator PO Box 8502 Madison WI 53708-8502 Phone: (920) 946-6909 E-mail: exam-coordinator@acbr.org CURRENT ACBR BOARD MEMBERS Tawnia Adams, DC, DACBR President E-mail: president@acbr.org Christopher Smoley, DC, DACBR Secretary E-mail: secretary@acbr.org Alisha Russ, DC, DACBR Member-at-Large E-mail: aruss@acbr .
