### Create a StudySoup account

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

Already have a StudySoup account? Login here

# ADVNUMERICALANALYSIS MATH693B

SDSU

GPA 3.64

### View Full Document

## 85

## 0

## Popular in Course

## Popular in Mathematics (M)

This 135 page Class Notes was uploaded by Miss Gladys Lubowitz on Tuesday October 20, 2015. The Class Notes belongs to MATH693B at San Diego State University taught by P.Blomgren in Fall. Since its upload, it has received 85 views. For similar materials see /class/225267/math693b-san-diego-state-university in Mathematics (M) at San Diego State University.

## Popular in Mathematics (M)

## Reviews for ADVNUMERICALANALYSIS

### 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: 10/20/15

Numerical Solutions to PDEs Lecture Notes 12 Systems of PDEs in Higher Dimensions 2D and 3D Time Split Schemes Peter Blomgren ltblomgrenpeter gmail com Department of Mathematics and Statistics Dynamlcal Systems roup Computational Sclences Raearch Center San Diego State University San Diego CA 921827720 httpterminussdsuedu Spring 2010 2D and JD Time Split Schemes 539 123 Outline 0 Recap 0 Last Time 9 Beyond 1D space 0 Mostly Old News with some Modifications 0 Instabilitites a Synthetic Example 0 Multistep Schemes a Finite Difference Schemes 0 The Leapfrog Scheme 0 The AbarbaneI Gottlieb Scheme 0 More General Stability Conditions 6 Time Split Schemes m 20 and 3D Time Split Schemes 7 223 Last Time The World is not One Dimensional Discussion Lower Order Terms and Stability Proof Dissipation and Smoothness Example Crank Nicolson in Dissipative Mode p fixed 0 o 0 Example Crank Nicolson in Non Dissipative Mode A fixed 0 0 Boundary Conditions accuracy ghost points 0 Convection Diffusion Grid restrictions due to the physics Reynolds or Peclet number of the problem upwinding 2D and JD Time Split Schemes m 323 In order to model interesting physical phenomena we often are forced to leave the confines of our one dimensional toy universequot The good news is that most of our knowledge from 1D carries over to 2D 3D and nD without change Such is the case for convergence consistency stability and order of accuracy The bad news is that the analysis necessarily becomes a little messier we have to Taylor expand in multiple dimensions all of which will affect stability etc From a practical standpoint things also get harder the computational complexity grows we go from On to 0 nd spatial grid points 5 More computations more storage more challenging to visualize in a meaningful way 539 2D and JD Time Split Schemes 7 423 Moving to Higher Dimensions Moving to Higher Dimensions Stability 1 of 2 We start out by discussion stability for systems of equations both hyperbolic and parabolic and then move on to a discussion of these systems in 2 and 3 space dimensions There is some news in testing for stability instead of a scalar amplification factor g0 we get an amplification matrix We obtain this matrix by making the substitution Wquot w Gquote 399 The VeCtOi VerSiOhS Of our medel Preblems are 0f the form The stability condition takes the form VT gt O 3C7 such that for 0 g nk g T we have a AUX 0 a B xx Hall CT where G is a d vector and the matrices A7 B are d gtlt d A must be Computing the G to the nth power may not be a lot of fun for a large diagonalizable with real eigenvalues and the eigenvalues of B must matrix G For hyperbolic systems this simplifies when G is a have Positive real Part polynomial or rational function of A this occurs in the Lax Wendroff and Crank Nicolson schemes There is very little news here for instance The Lax Wendroff In this case the matrix which diagonalizes A also diagonalizes G and scheme for the vector one way wave equation and the the stability only depends on the eigenvalues a of A eg for Crank Nicolson schemes for both vector equations look just as in Laxwendroff we must have Jaix g 1 for i 17 7d the 1D case but with the scalars a7 b replaced the matrices A7 B m m 2D and 3D Time Split Schemes 7 523 2D and 3D Time Split Schemes 7 623 Moving to Higher Dimensions Stability 2 of 2 Example An Unstable Scheme 1 of 2 For parabolic systems especially for dissipative schemes with p we COhSlder the someWhat artlflClal example constant similar simplifying methods exist U1 0 The matrix which transforms B to upper triangular form J LIZ J J O J B g lBU can also be used to transform G to upper triangular t form G Then if we can find a bound on HG a similar bound and the first order accurate scheme applies to n1 7 n 7 n 7 n n For more general schemes the situation is more complicated A V777 T Vquot 6Wm1 2W Wmil necessary condition for stability is Wm T Wm The corresponding amplification matrix is lgul S 1 Kk 2 6 for all eigenvalues g of G However this condition is not G J 1 4E 539 5 J sufficient in general 0 1 m 5353 2D and JD Time SpliL Schemes 7 723 2D and 3D Time Split Schemes 7 323 Example An Unstable Scheme 2 of 2 Multistep Schemes as Systems 1 of 2 The eigenvalues of G are both 1 but We can analyze multi step schemes by converting them into systems form eg the scheme 6 Gquot 1 4nesm2 E O 1 1 K vquot 5 Zaxevws Hence HG 7r2H On which shows that the scheme is Vo unstable D can be written in as a K 1 system The good news is that the straight forward extensions of stable schemes for single equations to systems usually results in stable VH1 GM 9 schemes As for scalar equations lower order terms resulting in Ok where V 17 5 Vn K T The matrix C67 is the modifications of the amplification matrix do not affect that companion matrix of the polynomial with coefficients 7a5 stability of the scheme given by 539 5351 2D and 3D Time Split Schemes 7 923 2D and 3D Time Split Schemes 7 1023 Multistep Schemes as Systems 2 of 2 Some Comments 30 31 H aKil 3K For scalar finite difference schemes the algorithm given in the 0 H 0 0 context of simple von Neumann polynomials and Schur C61 0 H 0 0 polynomials is usually much easier than trying to verify an estimate like llGnll CT 0 0 I 0 For multi step schemes applied to systems of equations there is no working extension of the theory of Schur polynomials so We note that this form of the companion matrix seems to be writing the scheme in the form of a one step scheme for an somewhat non standard both PlanetMath0rg and enlarged system is usually the best route in determining the mathworldwolframcom give a slightly different but essentially stability for such schemes equivalent form m 5393 2D and 3D Time Split Schemes 7 1123 2D and 3D Time Split Schemes 7 1223 Finite Difference Schemes in Two and Three Dimensions Leapfrogging Along in 2D 1 of 3 As stated earlier our definitions for convergence consistency and stability carry over to multiple dimensions however the von Neumann stability analysis becomes quite challenging We consider two examples First we consider the leapfrog scheme for the system tA X Buy0 where A B are d gtlt 039 matrices We write the scheme n1 7 n71 Vm VZm 2k A60X 560 ng 0 where 60X 60y are the central space differences in the X and In order to perform the stability analysis we introduce the Fourier transform solution V V 1 2 formally we can let vz jm w Gne wle mez where 9 h i 12 With A1 khl and A2 khg we get the recurrence relation VH1 2391Asin61 Angin62V i VH o ie we are interested in the amplification matrix G which satisfies G2 2 AlAsin61 Angin62 G i o The scheme can be rewritten as a one step scheme for a larger system and we can derive an expression for G for that system and check l Gnll 3 CT However it is very difficult to get reasonable y directions respectively 99 conditions without making some assumptions on A and B mm ID and 3D Time Split Schemes 7 1323 2D and 3D Time Split Schemes 7 1423 Leapfrogging Along in 2D 2 of 3 Leapfrogging Along in 2D 3 of 3 The most common assumption which rarely has any connection to 13 reality is that A and B are simultaneously diagonalizable That is we assume there exists a matrix P for which both PAP 1 and PBP 1 are diagonal matrices We let 04 and u be the A diagonal entries of these matrices and note that with the linear transform v39v P7 we get 039 uncoupled scalar relations W3 2 Ala sin61 A25 sin62 wquot 7 mil o m where 1 1 d This is somewhat more tractable we can The mOSt PeSSimiStiC Stability region is glVe by reuse our previous knowledge and we can conclude that the Allalmax A2l lmax lt 1 scheme IS stable if and only if where lozlmaX and l lmaX are computed from the separate AllaullAZl l lt 1 V1d m diagonalizations oannd B 5353 2D and 3D Time Split Schemes 7 1523 2D and 3D Time Split Schemes 7 1623 The Abarbanel Gottlieb Scheme 1 of 2 The Abarbanel Gottlieb Scheme 2 of 2 A resource saving modification to the leapfrog scheme which Slnce ObVlOUSly allows for larger time steps is given by Way sin61cos62 A2 ysin02cos61 maX1laul2l ul lsin91llcos02llsin92llcos91l vquot1 7 vquot 1 v v v v 7m 2k 7m A60X lmll 2 lmTl 560y 17m 2 Tl m O S max1lauli Azl ul sin261 c0529112 sin202 c0529212 max1lal Agl yl With the simultaneous diagonalizable assumption the stability condition is given by The two conditions M104 sin61 cos62 A23 sin62cos61l lt 1 Al au lt 17 MW lt 1 A sequence of inequalities can make some sense out of this are suf cient for stabmty and aso necessary 99 5351 2D and 3D Time Split Schemes 7 1723 2D and 3D Time Split Schemes 7 1323 More General Stability Conditions Time Split Schemes 1 of 3 It is possible to derive more general stability conditions without simultaneous diagonalization lfthe problem is hyperbolic easiest argued from the physics then the matrix function A51 352 is uniformly diagonalizable ie we can find a matrix P with uniformly bounded condition number so that P A51 552 8 1 05 is a diagonal matrix with real eigenvalues The stability condition becomes max max DA sin l A sin l lt1 193 W l 1 1 2 2m Sometimes this can be done with reasonable effort in other cases it is a big task 2D and JD Time Split Schemes 7 1923 Much of the work when it comes to devising practically useful schemes in higher dimensions is in the direction of dimension reduction ie reducing the problem to a sequence of lower dimensional problems Consider 6 6 A7 Bi O l 6xl l 6yl One way to simplify this is to let A act with twice the strength during half of the timestep with 5 turned off and then switch ie 6 ut2A7uO t0 t t0k2 8X ut2B uO t0k2 t t0k 0y 5353 2D and JD Time Split Schemes 7 2023 Time Split Schemes 2 of 3 Time Split Schemes 3 of 3 The analysis of time split schemes becomes quite interesting to say the least o If we use second order accurate difference schemes the overall scheme is second order accurate only if the order of the splitting is reversed on alternate time steps References For More Details 1 D Gottlieb Strang type Difference Schemes for Multidimen sional Problems SIAM Journal on Numerical Analysis 9 1972 pp 650 661 2 G Strang On the Construction and Comparison of Difference o Stability for split time schemes do not necessarily follow schemes 5AM Journa on Numericai Anaiysisy 5 1968 Pp from the stability of each of the steps Only in the case where 506 517 the amplification factors if being matrices commute is this true see 1 and 3 RJ LeVeque and J Oliger Numerical Methods Based on Ad 39 39 ditive Splittings for Hyperbolic Partial Differential Equations 9 Prescribing appropriate boundary conditions is a challenge Mathematics of Computation 40 1983 PP 46g4g7 see 539 m 20 and 3D Time Split Schemes 7 2123 20 and 3D Time Split Schemes 7 2223 Homework 3 Due 3192010 1234pm GMCS 587 Strikwerda 632 Theoretical Strikwerda 6310 Numerical Strikwerda 6314 Theoretical 5331 2D and JD Time Split Schemes 7 2323 Peter Blomgren blomgren peter gmail com Department of Mathematics and Statistics ynamicai Systems Group Computationai Sciences esearch Center San Diego State University San Diego CA 9218277720 httptei 39minnssisuadu Spring 2010 Peter Blom Outline 0 Recap 0 The ADI Method 9 SecondOrder Time Equations 0 Examples WaveEquation 0 What39s in a name Hyperbolic vs Parabolic a Finite Differences 0 Stability for SecondOrder Equations 0 Example CTCS for the Wave Equation 6 Example Order22 for the EulerBernoulli Equation 0 von Neumann Polynomials and Stability Peter Blomgren Second Order Equations e Differences 7 223 Recap The ADI Method Last Time The ADI Method The Alternating Direction Implicit ADI method allows us to solve primarily parabolic equations in multiple space dimension by slicing higher dimensional problems into onedimensional sub problems The slicing pushes the boundary of what size problem is computationally feasible A fully discretized ADI scheme based on a Crank Nicolson iteration for Ur Alu l Agu um i uyy is given by k k k k I 7 54th I 7 5AM Vr39Tl I 34th I 3AM Vquot There are several approaches to solving this including the Peaceman Rachford and D39Yakonov schemes Peter Blomgren Second Order Equatiol le Differences 7 323 SecollltLOrltler von n rm My Second Order Time Equations We now turn our attention to PDEs with second order time derivatives eg Hg 7 azuxx O The wave equation utt l bqu O The Euler Bernoulli beam equation Hg 7 czutm l bqu O The Rayleigh beam equation Most of our previously developed methods and theory can be applied to these equations with minor modifications Most prominently the definition of stability must take the second order derivative in time into account Peter Blomgren Second Order Equau 39 fferences 7 423 SecoilltLOrltler Time Equations Dir von Noumalm Pm hility The Second Order Wave Equation In one spacedimension the second order wave equation is given by Hg 7 azuxx 07 where a is a non negative real value the speed of propagation The initial value problem for this equation requires two initial conditions typically given as u07 X u0x7 ut0x u1x With exact solutions given by utx 2 2a u0xiatu0xat 1 Xatu1sds X fat Second Order EquaLio s y SecoilltLOrltler Time Equations rm V Di d The Second Order Wave Equation The exact solution shows that we have two characteristic speeds ia associated with the second order wave equation In the Fourier domain the solution is given by Dtw Ugwcosawt lw u4rwe lam l E All these expressions show that the general solution consists of two waves one moving to the right and one moving to the left Second Order Equ SecollltLOrltler von Neumann peiyn The Second Order Wave Equation The exact solution shows that we have two characteristic speeds ia associated with the second order wave equation In the Fourier domain the solution is given by 3t7w Ugwcosawt lw 3weia t we ia t All these expressions show that the general solution consists of two waves one moving to the right and one moving to the left Another way of seeing this is to split the differential operator J 82 2 82 J J 8 8 J J 8 8 J a u l 37 u 0 7 7 7 i 7 37 i 2 2 81 8X 81 8X 81 8X Peter Blomgren Second Order Equau 39 fferences 7 523 Example WaveEqualion s SecondrOrder Time Equations 39u Dilfel r m No Im Po m y Example 1 The Second Order WaveEquation We fix a 1 and use the following initial data cos 7rx 2 X lt17 u0x O lxlgL u1x0 The exact solution is given by unsegim whm w i Figure Snapshots at T 0 T 05 and T 10 See also the movie m waveOlmpg Second Order Equ SecondrOrder Time Equations 6 nit m Noulvialm Poly ity Examples WaveEqualion Example 2 The Second Order WaveEquation We fix a 1 and use the following initial data 1 le S 17 000707 1X O X gt17 The exact solution is given by Xt utx u1s d5 lengthx 7 tx 139 7171 X7 Figure Snapshots at T 0 T 05 and T 10 See also the movie m wav902mpg Second Order Equa SecoilltLOrltler Time qualions Exam 399 WZWE quotmoquot 39 Dinbi s m r m No m m y Example 3 The Second Order WaveEquation We fix a 1 and use the following initial data 7 7COS7TX2 lxl g 17 7 1 lxl g 17 lxl lo ixigt1 7 ll 0 ixigt1 The exact solution is given by 71 we utx 7 7 cos 7 fl cosltwlxgtgt lengthxi mmm 7111 Figure Snapshots at T 0 T 05 and T 10 See also the movie m waveO3mp Second Order Equ isrm Wlalslnaname 7 Hyperbolic vs Parabolic y lof2 SecomLOrder Time Equations Differ 7 von m P d Hyperbolic vs Parabolic The names hyperbolic and parabolic are historical and originate from the fact that the symbols of the second order equations are similar to the equations for hyperbolas and parabolas Laplace transforming in time and Fourier transforming in the spatial coordinates and setting 7 i5 gives For Hg 7 azuxx 0 the symbol is given by s2 7 32772 For Ht 7 buxx 0 the symbol is given by s 7 bnz The solutions 52 7 32772 C17 5 7 b72 C2 describe hyperbolas and parabolas respectively Second Order Equ What39s in a name 7 Hyperbolicquot vs Parabolic SecondrOrder Time Equations nite Ditfei s m No Im Po id Hyperbolic vs Parabolic 2 of 2 sZamhc1 sbm2c2 Figure The hyperbolas 27372 C1 in the 7157plane LEFT and the parabolas s 7 b72 C2 in the 7157plane RIGHT The names are just based on this formal similarity but are now fixtures in the language of PDEs The key to hyperbolic systems is that the solution propagates with finite speeds and the key to parabolic systems is that the solution becomes smoother than its initial data Second Order Equ vonNoI mun PLlyl 7 Hyperbolic V Parabolic The Euler Bernoulli Equation 1 of 2 The EulerBernoulli equation Ur bz llme describes the vertical motion ofa thin horizontal beam with small displacements from rest Using the Fourier transform it is straightforward to write down the exact solution 1 r A A sin bwzt utx 7 5quot u0 cosbw2t u1w72 dw x27r 700 bw iwxbwt w exibwt iou dw 1 oo e x27r 700 The second formulas shows that the propagation speed is ibw hence the equation is dispersive Second Order Equation Finite Difference SecondrOrder Time Equations V m Non mm My 3quotquot Whalsinaname 7 Hyperbolicquot vs Parabolic The Euler Bernoulli Equation 2 of 2 The Euler Bernoulli equation does not have finite speed of propagation bw is unbounded hence it is not hyperbolic Further there is no increased smoothness in the solution as time evolves hence it is not parabolic Second Order Equation Finite Difference Our definitions for convergence consistency and order of accuracy remain the same however we must modify or definition of stability minimum Kerr my warm l A finite difference scheme thv O for an equation that is secondorder in t is stable in a stability region if there is an integer J and for any positive time T there is a constant CT such that 00 J 00 h 2 lv2l2139m2CThZ Z w m7oo j0m7oo for all solutions v and for 0 g nk g T with k7 h E A The factor 1 n2 is new and shows that we allow a linear growth in t J is almost always 1 since data must be given at two timelevels Peter Blon ren bliiiiirf vift m e Band Qr ier Equations In the van Neumann anaiysis we must require that the at least two amplification factors satisfy igyi S 1 Kk If there are no Tower order terms then the stability condition is igyi g 1 with do we roots on the unit circle c1 TF1 WW to Wool Mr i If the ampli cation polynomial lt1gtg7 6 for a second order timedependent equation is explicitly independent ofh and k then the necessary an suf cient condition for the finite difference scheme to be stable is that all roots 6 satisfy the following conditions a igy9i S1 and b ifigyw 1 then igyw must be at most a double root imliih Maum TCS for the Wave Equation u in Example I rence Von No Ialln PLIyl mi Stability Example Central Time Central Space The standard second order accurate scheme for ut azuCY is V2 2V2 V271 2 V344 2V2 V271 k2 a h As usual we set v w gneime and factor out common terms to get g 7 2 7 g 1 7432A2 sin2 gl2 7 g lZY i2iaAsin 32 gl2 7 g l2 i2iasin O gi2l aAsing12 71 giZ iiaAsin i 17 32A sin2 Flnile Difference gi lt 1 7 azAz sin2 i iaAsin Second Order Equation Exam ple Centra l Time Central Space We have 2 t9 t9 gi lt 17 aZ xzsin2 i iasin and it is clear that as long as a g 1 we have lgil 1 At 9 0 gr g The equality also occurs when a 1 and t9 7139 Since we can allow two equal roots on the unit circle we have shown that the scheme is stable if and only if a g 1 Note Usually we take a lt 1 to avoid the linear growth of the wave with j 7139 where gi if2 71 even though this growth formally does not affect the stability of the scheme Peter Blomgren Second Order Equau 39 fferences 7 1723 rences Von Non ialm PLlyl ml Stability L A i i i w Example 722 for die EulerrBemoulli Equation Example The Euler Bernoulli Equation The simplest secondorder accurate scheme is given by n1 n n71 vm 7 2vm vm n n n n n ibz Vm2 4quoter1 6quotm Illmil Vm72 k2 4 The amplification factors are given by the roots of 6 k g72g71 716b2M2sin4 M7 It is stable if and only if 2szin2ltggt 1 gt ng Second Order Equation Finite Difference st rences m rug mun My ml Stability Example 422 for me EulerrBernoulli Equation Getting Started Computing v All schemes for second order time equations require some initialization of vi With u07X 0007 MOM U1X given the simplest procedure is based on the Taylor expansion 1 uk7 X N u0x ku1x EkzunOX 0 k3 Using u azum we get a second order accurate initialization from v1 7 u azk 2 7m kl 0 u11mT Xv The initialization should be of the same order of accuracy as the Note scheme in order not to degrade the overall method m Second Order Equation Finite Difference Secm lrdel39 lTimel El Finite nin ei Von Neumann Polynomials and Stability mam nidifStability We can modify our previously defined algorithms for von Neumann and Schur polynomials to test for the stability of secondorder schemes First we extend Id Def on Von Neumann Polynomial The polynomial LP is a von Neumann polynomial if all its roots rp satisfy lrpl S 1 Emma Cm Mam await de l The polynomial go is a von Neumann polynomial of order q if all its roots ry satisfy the following conditions a lryl g1 and b the roots with lryl 1 have multiplicity at most q A von Neumann polynomial of order 0 is defined to be a Schur polynomial PelerBlon ren i r P Secondr rdel39 lTimel El inite niirei Von Neumann Polynomials and Stability Old Theorem Von Neumann Polynom39al Te t pd is a von Ne mann polynomial of degree d if and only if either a ltpd0l lt and Lpdsl is a von Neumann polynomial of degree d7 1 or b Lpdsl is identically zero and Lp d is a von Neumann polynomial Wm CW Newman WWW Tree s l A polynomial god ofexact degree d is a von Neumann polynomial of order q ifand only if either a V lgod0l lt lgo0l and godsl is a von Neumann polynomial of degree d 7 1 and order q or b godsl is identically zero and gag is a von Neumann polynomial of order 7 Peter Blon ren Secom Finite Von Neumann Polynomials With this more general definition we have that a simple von Neumann polynomial is a von Neumann of degree Also these generalizations explain some of the parts of the algorithm 94me l Start with Lpdz of exact degree d and set NeumannOrder 0 while d gt 0 do 1 Construct LPZ 2 Define ad MAO 7 mam ltgt 3 Construct the polynomial 1M2 Lp04pdl 7 pd04pl 4 1 If 1M2 E 0 then increase NeumannOrder by 1 and set Lpd1z 4412 4 2 Otherwise if the coefficient of degree d 7 1 in 1M2 is 0 then he polynomial is am Neumann polynomial of any order terminate algi 4 3 Otherwise set Lpd1z 1M2 endwhile decrease d by 1 4 Enforce appropriate Conditions on 2 seconrurraer i 39 Von Neumann Polynom von Neumann Polynomials and Stability The new theorem and the algorithm can be used to analyze stability for second order equations If Clgtg0 is the amplification polynomial of finite difference scheme for a second order equation for which the restricted condition lgl0l 1 can be used the the scheme is stable if and only if Clgtg0 is a von Neumann polynomial of order 2 Next time Boundary conditions two and three spatial dimensions Peter Blomgren Second Order Equau 39 fferences 7 2323 Numerical Solutions to PDEs Lecture Notes 9 Dissipation and Dispersion Peter Blomgren blomgren peter gmail com Department of Mathematics and Statistics Dynamicai Systems Group Computationai Sciences Research Center San Diego State University San Diego CA 9218277720 httpterminussdsuedu Spring 2010 Peter Blon ren Outline 0 Recap 0 Stability and von Neumann Polynomials O Hyperbolic PDEs Consistency Accuracy Stability a Dissipation 0 Introduction Leapfrog 84 LaxWendrofF 0 Adding Dissipation a Dispersion 0 Introduction 0 Dispersion Relations Peter Blomgren 53 mgvl u ly We developed a framework and in the end an algorithm for checking the stability of a general multistep scheme Such a scheme is stable exactly when its amplification polynomial 4p is a simple von Neumann polynomial Simplevom Neumann Pawnquot l The polynomial 4p is a simple von Neumann polynomial if 4p is a von Neumann polynomial and its roots on the unit circle are simple roots The derived algorithm involves comparing the magnitude of the first and last coefficients of the polynomial and then forming a polynomial of on degree less and recursively applying the same test to this polynomial Peter Blomgren 7 Recap u LI 1 i i h quotl in a Hyperbolic PDEs Conslslency Accuracy Stability Looking Back Our discussion of finite difference schemes for hyperbolic PDEs is almost complete Thanks to the Lax Richtmyer equivalence theorem we only need to consider schemes that are both consistent and stable From the discussion last time we have a very general framework for checking stability and we get consistency as a sideeffect of finding the order of accuracy of the scheme Usually finding the order of accuracy comes down to a Taylor expansion and in the hardest cases we can fall back and use the symbols of the scheme and PDE and congruence to zero for homogeneous equations in order to mechanize the analysis For derivation of high order accurate schemes we have a symbolic cacuus with the difference operators 6 6 60 62 etc m Peter Blomgren Di patioquot ion D Dissipation and Dispersion We now look at two final topics in the context of hyperbolic equations Dissipation and Dispersion Dissipative schemes Damping out high frequency waves which make the solution too oscillatory Recall the problem of parasitic modes of the Leapfrog scheme Dispersion Numerical Dispersion Refers to the fact that finite difference schemes propagate dif ferent frequencies at different speeds This causes the solu tion to change shape spread out as t grows Peter Blomgren alion and Di Introduch Leapfrog amp LaerendrofI39 Dissipation The order 22 Leapfrog scheme performs better than the order 12 Lax Friedrichs scheme but the solution tends to contain small oscillations Introduclio Leapfrog amp LaerendrotT The Lea pfrog Scheme Consider the leapfrog schemes with oscillatory initial conditions v 1 7n7 j 01 lnl ltlt1 It is straight forward to see that the leapfrog solution turns out to be VPquot 1 n Figure The leapfrog scheme propai gates the oscillatory solution This means that if we ever due to numerical error bad boundary conditions round off error etc introduce oscillations they never go away Introduch Leapfrog amp LaerendrotT L 1 n i i The Lax Wendroff Scheme Next we consider the Lax Wendroff scheme with the same oscillatory initial conditions v2 71m 7 71a Lax Wendrorrzma T20 Figure Here we can see that the magnitude of the highifrequency oscillations have m been dampened down to N 10 16 at time 239 Introduction Leapfrog 8 Laerendrofl39 The LaX Wendroff Scheme 2of3 The exact solution of the Lax Wendroff scheme with this particular initial data is given by v 17 232A2quot71mquot and since lal lt 1 here 08 we get rapid exponential decay of the oscillations The Lax Wendroff scheme is dissipative Peter Blomgren alion and Di Introduction Leapfrog x LaerendrofI39 The Lax Wendroff Scheme Definition Dissipative Scheme A scheme is dissipative of order 2r if there exists a positive constant c independent of h and k such that each amplification factor 59 satisfies lgy6l g 17 csinz lt6 9 2 lt 7 Zr 7 2 ltgt lgy6l 71 CSln lt2 The amplification polynomial for the LaxWendrofF scheme satisfies and is dissipative of order 4 for O lt laAl lt 1 Peter Blomgren alion and Di Di patioquot D ion Tm lQXthsm li rdaf A scheme is dissipative of order 2r if there exists a positive constant c independent of h and k such that each amplification factor 59 satisfies 6 lgy6l g 17 csinz ltgt lgy6l2 g 17 c sinz The amplification polynomial for the LaxWendroff scheme satisfies 6 lg 9l2 17 432A217 aA sin4 iue of order 4 for O lt laAl lt 1 The leapfrog scheme and the CrankNicolson scheme are strictly mam dissipative since their amplification factors are identically 1 in magnitude The LaxFriedrichs scheme is nonstrictly nondissipative Peter Blomgren Di pal D n io io Adding Dissipation Adding Dissipation to Schemes It is possible to add dissipation to non dissipative schemes however care must be taken so that we do not affect to order of accuracy of the scheme The modified Leapfrog scheme given by n71 fn vrrr39fl 7 vrquot1 vrrr391 7 vr il e h6 4 vm m 2k 3 2h 7 is second order accurate for small values of e the corresponding amplification factor is gi0 fiasint9 i 1 32V Sinzw 7 ESin4 alion and Di Peter Blomgren Adding Dissipation to Schemes For small enough 6 the magnitude of the roots are given by 0 lgi0l2 17 esin4 hence the modified leapfrog scheme is dissipative of order 4 Figure Comparison of LEFT modified leapfrog 5 05 and RIGHT standard m leapfrog A 05 h 0V1 Li H Introduction Di 0 Dispersion Dispersion Introduction We write the solution to the homogeneous one way wave equation using the Fourier inversion formula 1 00 u tx 7 e quote at d 7 mm L 3tw from this representation we can identify t 19w e mk 31 w DiSSiI39S39ILiOD Ispevsuon Dispersion Introduction We write the solution to the homogeneous one way wave equation using the Fourier inversion formula 1 00 u tx 7 e quote at d 7 mm L 3tw from this representation we can identify t 19w e mk 31 w If we consider a one step finite difference scheme the propagation in Fourier space is given by W1 gumquot Peter Blomgren DlSSlI39S39ILlOD Ispevsuon Dispersion Introduction We write the solution to the homogeneous one way wave equation using the Fourier inversion formula 1 00 u tx 7 e quote at d 7 l milm L 3tw from this representation we can identify t 19w e mk 31 w If we consider a one step finite difference scheme the propagation in Fourier space is given by W1 gumquot Clearly if the scheme is good then gh N e igak m Peter Blomgren Titian Introduction his r Di Dispersion Dispersion Introduction In order to clearly show the connection between the scheme gh and the solution of the PDE e iwak we write gm igh i WWW where 04h5 is interpreted as the phase speed the speed at which waves of frequency 5 are propagated Peter Blomgren alion and Di DlSSlI39S39ILlOD Ispevsuon Dispersion Introduction In order to clearly show the connection between the scheme gh and the solution of the PDE e quot we write gh lgh l WWW where 04h5 is interpreted as the phase speed the speed at which waves of frequency 5 are propagated lf 04h5 E a then all waves would propagate with correct speed For many finite difference schemes this does not hold and the difference a 7 04h5 is known as the phase error Peter Blomgren 39 alion and Dispersion DlSSlI39S39ILlOD Ispevsuon Dispersion Introduction In order to clearly show the connection between the scheme gh and the solution of the PDE e quot we write gh lgh l WWW where 04h5 is interpreted as the phase speed the speed at which waves of frequency 5 are propagated lf 04h5 E a then all waves would propagate with correct speed For many finite difference schemes this does not hold and the difference a 7 04h5 is known as the phase error The phenomenon of waves of different speeds traveling with different speeds is known as dispersion Peter Blomgren 39 alion and Dispersion DlSSlI39S39ILlOD Ispevsuon Dispersion Introduction In order to clearly show the connection between the scheme gh and the solution of the PDE e quot we write gh lgh l WWW where 04h5 is interpreted as the phase speed the speed at which waves of frequency 5 are propagated lf 04h5 E a then all waves would propagate with correct speed For many finite difference schemes this does not hold and the difference a 7 04h5 is known as the phase error The phenomenon of waves of different speeds traveling with different speeds is known as dispersion The effect of dispersion is shapedistortion of the solution Peter Blomgren 39 alion and Dispersion n Introduction 0 Dispersion Dispersion Shape Distortion Lax Fried richs 208 T20 o dx110 dx1l20 39ILion Di i i Dispersion Dispersion Relations Dispersion Relations With 9 hE we write g0 lg0l e mu WW and dusting off some complex analysis we identify lm lg0l tan laiowl m39 When lg0l 1 we have sin a00 lm g0 Peter Blomgrei Dissipation and Dispersion 39ILion Di i i Dispersion Dispersion Relations Example The Lax Wendroff Scheme The amplification factor Tor the Lax Wendroff scheme is g0 17 2a2sin2 7 iasin0 By the preceding arguments we have tan a00 With a little help from our friend Taylor this can give us some information about 30 5M9 9 1 7 92 0 94 tan9 9 1 92 0 94 tanT19 9 1 7 9 0 94 m Peter Blomgrei Dissipation and Dispersion DlSSlI39S39ILlOD ispeision Dispersion Relations Example The Lax Wendroff Scheme We get 1 aw a 17 602 17 a2 0 04 0 When 9 is small and lal lt 1 then 040 lt a the numerical solution will tend to trail the exact solution When lal z 1 then the dispersion phase error is smaller Peter Blomgren 39 alion and Dispersion 39ILion Di i i Dispersion Dispersion Relations Example The Lax Wendroff Scheme For larger values of 0 the Taylor expansions do not hold and we must consider the full expressions g0 17 2a2sin2 7 iasin0 asin0 t 9 A0 an la l 172a2sin2 g When 9 7139 g 17 232 o If lal gt then g lt 0 and a7r 1 o If lal lt then g gt 0 and a7r 0 By consistency 040 a Peter Blomgrei Dissipation and Dispersion mien v Dispersion Dispe on Relations Example The Lax Wendroff Scheme Figure The phase speed for the LaX7WendrofT scheme 7 LEFT 069 lt AND RIGHT A 099 080 071 gt i a 1 DlSSlI39S39ILlOD Dispeision Dispersion Relations Rules of Thumb o For hyperbolic PDEs is it best to take lal as close to the sta bility limit as possible this keeps the dissipation and dispersion small eg the phase speed for Lax Wendroff for a 1 and 099 If we are interested in a propagating a particular frequency 5 then we must choose h so that h lt 7r Think Nyquist sampling m 39 alion and Dispersion Peter Blomgren Rules of Thumb Di L on Dispersion i Disper n Relations When the leapfrog and Lax Wendroff schemes are applied to the homogeneous one way wave equation with lal 1 there is no dispersion error These are exceptional special cases when atx is variable or the system is non trivial in any other way there is dispersion Dissirntion i Ispelslon Disperslon Relations Rules of Thumb c When the leapfrog and Lax Wendroff schemes are applied to the homogeneous one way wave equation with lal 1 there is no dispersion error These are exceptional special cases when atx is variable or the system is non trivial in any other way there is dispersion o The phase error is always an even function of 0 hence if a scheme has order of accuracy p then the phase error is of order p if p is even and order p1 if p is odd Peter Blomgren 1 i Disperslon Relations Rules of Thumb When the leapfrog and Lax Wendroff schemes are applied to the homogeneous one way wave equation with lal 1 there is no dispersion error These are exceptional special cases when atx is variable or the system is non trivial in any other way there is dispersion The phase error is always an even function of 0 hence if a scheme has order of accuracy p then the phase error is of order p if p is even and order p1 if p is odd When choosing a scheme for a particular application the amount of dissipation and dispersion can and should be used to choose between schemes see eg D Durran The Third Order Adams Bashforth Method An Attractive Alternative to Leapfrog Time Differencing Monthly Weather Review 119 1991 pp 702 720 Peter Blomgren 39 alion and Dispersion 5351 39ILion Di 1 i Dispersion Dispersion Relations Wrapping Up and Looking Forward That ends our discussion of finite difference schemes for hyperbolic PDEs Next time we start looking at parabolic PDEs of which the onedimensional heat equation 8 82 EUU 7X 7 btx ut7 X uioix MX is the simplest example Peter Blomgren alion and Di Nyquist Sampling The following result is fundamentally important in signal processing information theory WWW Sampling i In order for a band limited ie one With a zero power spectrum for frequencies 1 gt B baseband 1 gt 0 signal to be reconstructed fully it must be sampled at a rate 1 2 2B A signal sampled at 1 2B is said to be M qufs t 5521 9127 and 1 2B is called the L ency No information is lost ifa signal is sampled at the Nyquist frequency and no additional information is gained by sampling faster than this rate An audio CD is sampled at 44100 samplessecond allowing for reconstruction of signals up to 22050 Hz There is no evidence that human beings are sensitive to audio frequencies above 20 kHz and most people over the age of 35 are unable to hear sounds above 15716 kHz at 72 dB Peter Blon ren Numerical Solutions to PDEs Lecture Notes 10 Parabolic PDEs Overview Exact Solutions and Boundary Conditions Peter Blomgren ltblomgrenpeter gmail comgt Department of Mathematics and Statistics D namlcalSystems roup Computational Sclences Raearch Center San Diego State University San Diego CA 921827720 httpterminussdsuedu Outline o Parabolic PDEs 0 Introduction Model the 1D Heat Equation 0 Systems and Boundary Conditions 6 Finite Difference Schemes 0 Mostly Old News Spring 2010 539 m Parabolic PDEs 7 127 Parabolic PDEs 7 227 Parabolic PDEs The 1D Heat Equation 1 of 6 Our sim lest model for arabolic PDEs is the one dimensional heat equationp P We Fourier transform in the spatial coordinates the 1D heat equation with constant b and get a friendly ODE Ll b Ll XX 3 7bw23 mob 30w u0X u0X where the heat conductivity b must be non negative In the The SOlUtlon m the Fourler domam 395 glven by Simplest case b Is a constant but we can allow bt X 2 0 mtyw eibwzt 300d Parabolic PDEs show up in modeling of gas and fluid flow economic modeling the Black Scholes equation with a negative and by the Fourier inversion formula we get b and diffusion processes think pharmaceuticals spreading in 1 00 the body epidemiological studies animal of the year flu etc ut7X 7 e wxe bwzt ow dw x27r As with the one way wave equation in the hyperbolic case we can 00 learn a lot from studying the 1D heat equation and numerical in the physical domain solutions thereof for the parabolic case 535 Parabolic PDEs 7 327 Parabolic PDEs 7 427 The 1D Heat Equation 2 of 6 The 1D Heat Equation 3 of 6 The factor e bwzt in the previous two expressions shows that Llt7 X is obtained from Ll0X by damping the high frequency oo 00 content of u ut X eiwxiyeibw2t dw 0y any 039 Mr 700 w 700 Hence the solution operator for a parabolic equation Is a 00 1 dISSIpative operator 7 4 bt e X Y24bt 00 dy By using the Fourier transform formula W TOO 1 00 This expresses Llt7 X as a weighted average of Llo with a Gaussian A 7iw HOW 9 y 0y dy weight function 7139 700 e7x7y24bt G trxr y In the preVIous expreSSIon for Llt7 X we get 47rbt 1 00 I W2 1 00 I It has the property that Vt 2 0 LI tyx 7 e wxe t ei wyu d dw x27r 00 x27r 700 00 y 1 7x7y24bt e dy 1 V 47rbt 700 Next we Interchange the order of Integration dy lt gt dw m m Parabolic PDEs 7 527 Parabolic PDEs 7 627 The 1D Heat Equation 4 of 6 The 1D Heat Equation 5 of 6 T m I 1 Given the representation um 71 00 W4 uoy dy 47rbt foo Hiram we make the following observations 0 Llt7X is infinitely differentiable in t and X for t gt 0 Figure The Gaussian weight func 0X 2 039 then t7X 2 039 tion at time t 00017 0 017 01 y 0 b 1 are fixed 5331 5393 Parabolic PDEs 7 727 Parabolic PDEs 7 327 The 1D Heat Equation 1 O 1 6of6 Figure The exact solutions for the heat equation with initial condition 1 for Parabolic Systems and Boundary Conditions A system of the form at B xx AUX C1quot Ftx where U is a vector and A B and C are matrices is parabolic if RealB gt O B does not have to be positive definite nor symmetric nor must the eigenvalues be real there are no restrictions on A and C lxl g 1 and 0 otherwise We clearly see how the heat spreads out this corresponds quite well with our physical intuition m m Parabolic PDEs 7 927 Parabolic PDEs 7 1027 Parabolic lVPs and Well Posedness 1 of 3 Parabolic lVPs and Well Posedness 2 of 3 Theorem Well Posed Parabolic lVPs Recall The initial value problem for the system Definition Well Posed lVP t B xx AUX CG Ft7 X The initial value problem for the first order partial differential I I equation Pu O is we Posed if for any time T 2 0 there exists a Is well posed39and actually a stronger estmate holds For each constant C such that any solution LltX satisfies T gt 0 there 5 a COnStant CT SUCh that 00 oo 00 t 00 lutXl2 dX CT luoxl2 dX l tXl2 dx l xsXl2 dxds 700 700 700 0 700 00 t 00 for0 t T CT llJOXl2ClX lF5Xl2dXClS 700 0 700 for O S t S T nu 3053 Parabolic PDEs 7 1127 Parabolic PDEs 7 1227 Parabolic lVPs and Well Posedness 3 of 3 Boundary Conditions for Parabolic Systems 1 of 2 The bound in the theorem is stronger than the one in the definition since it gives a bound on the derivative of u with respect to X in addition to a bound on u The bound on uX implies that the solution to the system is infinitely differentiable for t gt O The proof of the theorem is quite straight forward using the Fourier transform Parseval39s Theorem 3 7wZBiwA C 3t7w 87u2BiuACt ow l trwl Kea b 2 l owl l twl2dw m l owl2dw A parabolic system with 039 equations defined on a finite interval requires 039 boundary conditions at each boundary common forms include T0Ut75 b0 T1LIXt T2ut5 b1 where 5 E 39 thus specifying the temperature or the relation between the flux and temperature at the boundary Note TO 6 Rdoxd T1 T2 6 R0140 Eg a perfectly isolated refrigerator would have t 00 co letf O w2l3swl2 dw d5 Kg l owl2dw 0 m m Parabolic PDE 7 1327 Parabolic PDE 7 1427 Boundary Conditions for Parabolic Systems 2 of 2 Finite Difference Schemes for Parabolic Equations Our previous definitions given in the context of finite difference h f h b l39 t39 f 39 t Boundary conditions are said to be well posed if the solution of 5 lFteS 0rd yper 0 IC equa lons cl convirgtiniet hconzls elncy the PDE depends continuously on the boundary data expressed in 5 fl l l y an alcfzurzicy wereg nera enOUg a ey pp y Without modification to finite difference schemes for terms of the matrices T0 T1 and B we must require that the parabolic equations 039 gtlt 039 matrix We start by considering the forward time central space scheme for T 7 T0 the heat equation mil2 vquot1 7 v v 72v v is non smgular bw When To ldxd we have Dirichlet boundary conditions n n ma speci ed temperatures and when T1 ldxd and T2 O we have We get the amplification factor usmg our old trick vm w g e Neumann boundar conditions 5 ecified tem erature fluxes y P P g67libeI072eile k T h2 5331 5353 Parabolic PDE 7 1527 Parabolic PDE 7 1627 Stability of the Forward Time Central Space Scheme With p we get g6l 1 7 41mm2 1 so that g6l g 1 as long as by g 5 From the form of g6 we immediately see that the scheme is dissipative of order 2 as long as O lt by lt Dissipation is desirable for parabolic equations since this implies that the numerical solution will become smoother as 139 increases mimicking the behavior of the PDE The stability condition expressed in terms p is generic for the parabolic universe of finite difference schemes p plays the role of A for hyperbolic problems 539 Some Forward Time Central Space Solutions On the next slides we see solutions to b 1 th LIXX X 6 733 139 E 0 T i 1 Xl S 1 07 o M gt 1 We have h g and k 0005 so that p 032 lt 05 The solutions are shown for 3 different time intervals 02 08 and 016 On slide 19 the solution is computed with boundary conditions Llt i3 0 On slide 20 the solution is computed with boundary conditions which match the exact solution this is quite artificial but allows us to compare the solutions see slide 21 m Parabolic PDEs 7 1727 Parabolic PDEs 7 1327 FT CS Solutions to the Heat Equation 1 of 3 FT CS Solutions to the Heat Equation 2 of 3 a it s l 7 r 479 7 cc 05 its a a 07 in 5 06 5 05 c 115 4 05 3 u a 04 0 as 2 2 u2 112 01 l I 3 2 1 u 1 2 a 9a 2 I n 1 2 3 Figure Solutions to the heat equation Flgure comparlscm f the lining for ranges T E 0 2 T E 0 8 and cal and exact solution in the lnflnlte l 7 domainH case shown for T 2 T 39 8 and T 16 m 5353 Parabolic PDEs 7 1927 Parabolic PDEs 7 2027 FT CS Solutions to the Heat Equation 3 of 3 More Schemes Backward Time Central Space The backward time central space scheme applied to th quX f is 1 1 1 7 vs b 7 W for hz m the amplification factor is 1 EM 29 Figure Comparison of the computed 1 4b 5 E and exact solutions to the heat equation 6 This implicit scheme is order 12 and is unconditionally stable When p 2 c gt 0 it is also dissipative of order 2 m m Parabolic PDEs 7 2227 Parabolic PDEs 7 2127 More Schemes Crank Nicolson The Crank Nicolson scheme applied to th quX f is given by 1 n1 7 n1 n1 n 7 n n Vi Vrquotn E Vm1 2quotm 1 Vm71 Vm1 2quotm Vm71 2 h2 h2 imi The corresponding amplification factor is 1 7 2bpsin2 g EM 29 1 2bp Sln 5 The Crank Nicolson scheme is implicit unconditionally stable and order 22 accurate When p is constant it is dissipative of order 2 m Parabolic PDEs 7 2327 More Schemes The Leapfrog Scheme The leapfrog scheme applied to th quX f is given by m b 2k The corresponding amplification polynomial is n1 n71 n 7 n n Vm T V Vm1 2Vm T Vm71 n 2 fin g62 8g6bpsin2 7 1 with roots we i 1 4bpsin2 4bpsin2 We see that lg6l gt 1 for most values of 6 hence the leapfrog scheme is unconditionally unstable Parabolic PDEs 7 2427 More Schemes The Du Fort Frankel Scheme 1 of 2 More Schemes The Du Fort Frankel Scheme 2 of 2 The Du Fort Frankel scheme is essentially a quotfixedquot leapfrog When 1 7 4b2p2 sin2 0 2 0 we get scheme 2bpl cos0l l 1 7 4b2p2 sin2 9 2b 1 n1 7 n71 n 7 n1 quot71 n lt lt M Vm Vm bvm1 Vm Vm Vm71 fn lgi0l7 l2bp 7 l2bp 2k h2 m 7 2 2 2 The corresponding amplification polynomial is and When 1 4b ll 539 0 lt 0 we get 2 2 2 2 2 2 2 1 2bpg6l 7 4bpcos6lg 7 1 7 2bp 57in S 2bpcosl9 4b p sm 0 7 1 4b p 7 1 lt 1 12bp2 4b2p2 4bp 1 with roots Thus this explicit scheme is unconditionally stable The one caveat is that it is onl consistent if k h tends to zero as h and k 0 to zero 2bp cos6 i 1 7 4b2p2 sin2 6 y g giw The caveat is a generic property for finite difference schemes applied to 1 219M parabolic problems as expressed in the theorem on the following slide 99 m Parabolic PDEs 7 2527 Parabolic PDEs 7 2627 Convergence of Explicit Schemes Theorem Convergence of Explicit Schemes An explicit consistent scheme for the parabolic system at B xx AUX CG Ftx is convergent only if k h tends to zero as k and h tend to zero Recall for hyperbolic problems Theorem There are no explicit unconditionally stable consistent finite difference schemes for hyperbolic systems of partial differential equations 5337 Parabolic PDEs 7 2727 Outline 0 Recap 0 Parabolic PDEs Stability Boundary Conditions 0 Schemes ForwardBackward Time Central Space Convection Diffusion Variable Coefficients SCl emeS3 crank39NlCOlsoni DU39FOVt Frankel a Stability Lower Order Terms 0 One step Schemes Numerical Solutions to PDEs Lecture Notes 11 Parabolic PDEs Peter Blomgren a Dissipation and Smoothness ltblomgrenpeter gma11comgt O Crank Nicolson Department of Mathematics and Statistics 6 Boundary Conditions Dynamical Systems Group Computational Science Raearch Center 0 2nd order one39SIdedi GhOSt Pomts San Diego State University San Diego CA 921827720 a Convection Diffusion o Numerlcs httptermlnussdsuedu 0 Upwi nd Differences Spring 2010 6 Variable Coefficients Stability 7 134 Stability 7 234 Last Time 1 of 3 Last Time 2 of 3 A quick introduction to parabolic PDEs Our model equation is Numerical SChemes for Ur bLIxx f the one dimensional heat equation Forward Time Central Space Exact solutions to the 1D heat equation in infinite space usmg the Fourier transform v3 7 VI 7 b V371 7 2V 1 vigil i 2 The solution corresponds to a damping of the high frequency k h content of the initial condition 5 the parabolic solution operator is dissipative fm Explicit stable when by S where p I772 order 12 dissipa tive of order 2 For t gt O the solution of the heat equation is infinitely differentiable Backward Time Central Space Since parabolic PDEs do not have any characteristics we need V37 V37 2Viilq1 n1 boundary conditions at every boundary Typically we specify Ll k 2 m fixed temperature Dirichletquot the normal derivative LIX Implicit unconditionally stable order 12 dissipative of order temperature flux Neumannquot or a mixture thereof 5351 Stability 7 334 Stability 7 434 Last Time 3 of 3 Lower Order Terms and Stability 1 of 2 CrankNicolson For hyperbolic equations we have the following result VH1 7 Vquot b Vn1 7 2Vn1 Vn1 vquot 7 2quot vquot 1 Theorem Stability of One Step Schemes m m 7 m1 m m71 m1 m m71 7 fquot1fmn I I k 2 h2 2 A consistent one step scheme for the equation Implicit unconditionally stable order 22 dissipative of order ut aux by 0 2 when p is constant is stable if and only if it is stable for this equation when Du Fort Frankel xed leapfrog b 0 Moreover when k Ah and A is a constant the stability v37 7 vl j fl v371 7 vfffl v fl 7771 n condition on gh 7 k h is 2k b h2 T fr lt Explicit unconditionally stable order 22 dissipative of order 2 lgwlolo 1 when p is constant It is only consistent if kh tends to zero as h and k go to zero Similar results do not always apply directly to parabolic m equations 5351 Stability 7 534 Stability 7 634 Lower Order Terms and Stability 2 of 2 Dissipation and Smoothness The fact that a dissipative onestep scheme for a parabolic equation generates a numerical solution with increased smoothness as t provided that p is constant is a key result so lets show that it is indeed The problem is that the contribution to the amplification factor from the first derivative is sometimes often 9 which is greater than 9k as k 0 For instance the forwardtime centralspace scheme applied to true at buxx 7 aux cu gives the amplification factor We Start with the following theorem 6 g 1 7 4busin2 7 iaAsin6 ck Theorem A one step scheme consistent with The term ck does not affect stability but the term containing A xkya cannot be dropped when u is fixed In this particular case we get u buxx 6 2 g z lt1 7 4b sin2 7 azk sinzw that is dISSIpatlve oforder 2 With it constant satisfies n n1 2 Z V 2 lt 0 2 and since the first derivative term gives an O k contribution to lglz it does not llV llz l Ck ll6v llz llV llz affect stability Strikwerda p149 This is also true for the backwardtime V1 centralspace and CrankNicolson schemes m for all initial data v0 and n 2 O 5 Stability 7 734 Stability 7 334 Dissipation and Smoothness Proof 1 of 2 Dissipation and Smoothness Proof 2 of 2 Proof Let 0 be such that lg19l2 g 1 7 0 sin2 dissipative scheme of order 2 Then by we have VVWE g0VVEi WNW g0l2lvV l2 WW 7 co WW By summing this inequality for 1 O7 n we get using 1 kh z 2 We get 1 lvn1lsll2 ck lf6V l2 ivolsllz 0 1 Integration over 5 and using Parseval s relation we get A c k quot 1 0 A A lvquot l2 07 gsm 5 WE lv l2 llvquot1ll ck ll6v ll llvOlli M VZO 110 NeXt we use which is the inequality in the theorem D 25in Q em 7 1 39 hlz h V V6VVE Here J denotes the Fourier transform m 5331 Stability 7 934 Stability 7 1034 Dissipation and Smoothness 1 of 2 Dissipation and Smoothness 2 of 2 We can use the theorem to show that solutions become smoother with time lt2 norms of the high order differences approximating high order derivatives tend to zero at a faster rate than the norm of LI Since lg6l S 1 we have llvquot1ll2 S llVVllg We note that 61v being a finite difference is also a solution to the scheme therefore we have H61 vquot1ll2 S ll6vquotll2 That is both the solution and its differences decrease in norm as time increases We apply the theorem and get llv 1ll ctll6v ll S llvolli which shows for nk t gt O that ll6vnll2 is bounded and we must have C ll6v ll 7llvoll 0 Stability 5351 7 1134 The argument can be applied recursively since 6 v satisfies the difference equations we find that for nk t gt O and any positive integer r that div is also bounded Thus the solution of the difference scheme becomes smoother as 139 increases It can be shown that if v 7 LltXm with order of accuracy p then div 7 6utnxm with order of accuracy p These results hold if and only if the scheme is dissipative 5353 Stability 7 1234 Example Crank Nicolson dX 120 dt 120 M 20 Example Crank Nicolson dX 140 dt 140 M 4O Figure The CrankNicolson scheme applied to the initial condition in panel 1 with zero Figure The CrankNicolson scheme here we have cut both h and k in half compared with flux boundary conditions We know that CrankNicolson is nondissipative if A remam the previous slide On the next slide we show the result of keeping it kh2 constantm constant see next slide which case the scheme is dissipative Stability 7 1334 Stability 7 1434 Example Crank Nicolson dX 140 dt 180 M 20 Example Crank Nicolson dX 180 dt 180 M 8O Figure The CrankNicolson scheme here we finally get some damping in the oscillations of the solution Dissipation is a convergence resulti Figure Surprlslneg reflnlnlng in X brings back the overshoot artefacts m Stability 7 1534 Stability 7 1634 Example Crank Nicolson dX 120 dt 140 M 10 Example Crank Nicolson dX 120 dt 180 M 5 Figure Coarsening in X dx 120 instead of dx 140 lessens the carrying capaci of highfrequency content of the grid h Figure m Stability 7 1734 Stability 7 1334 Boundary Conditions Again More Accurate Boundary Conditions 1 of 2 Since parabolic problems require boundary conditions at every boundary there is less need for numerical boundary conditions compared with hyperbolic problems We briefly discuss implementation of the physical boundary conditions Implementing the Dirichlet specified values at the boundary points boundary conditions is straight forward The Neumann specified fluxderivative is more of a problem for instance one sided differences 3utnxo N v1 7 v67 3utnXM N VA 7 VIC1 3X N h 7 3X N can be used but these are however only first order accurate and will degrade the accuracy of higher order schemes Boundary Conditions 5351 7 1934 Second order one sided accurate boundary conditions are given by ammo N 7v2quot4v1quot asvoquot aumxn N vim 74ml 3vii 3X N 2h 7 3X N 2h It is sometimes useful to use second order central differences and introduce ghost points for the boundary conditions eg a quottn7X0 N V1 7 V21 3X N 2h How is this useful Consider a given flux condition uxtXo Mt then V1 iv l 7 n 7 n 7 7 on lt2 v71 7 v1 72h n 2 7 5353 Boundary Conditions 7 2034 More Accurate Boundary Conditions 2 of 2 Now if we are leapfrogging Du Fort Frankel style the scheme can be applied at the boundary m O n1 n71 n7 n1 n71 n V0 V0 bVl V0 V0 V71 fr 2k h2 l m vs 7 b v1quot 7 ch V321 7 V17 2w 2k h2 m Ideas like these are commonly used 539 Boundary Conditions 7 2134 The Convection Diffusion Equation Numerics 1 of 3 First we consider the forward time central space scheme n1 7 n Vm V ma n 7 n n 7 n n Vm1 Vm71 i b Vm1 2Vm Vm71 h2 which is first order in time and second order in space Since stability requires by S 12 we must have k N h2 so the scheme is second order overall For convenience lets assume a gt 0 define p and 7 E 7 2A 39 a 7 2b 7 72W we can write the scheme as vs 7 1 7 21mm bM1 7 av1 bM1 7 com Based on previous discussion of parabolic PDEs we know that llut S Hut 7 if t gt t the peak value is non increasing m Convection7Diffusion 7 2334 The Convection Diffusion Equation Many physical processes are not described by convection transport eg the one way wave equation or diffusion eg the heat equation alone An oil spill in the ocean or a river is spreading by diffusion while being transported by currents the same goes for your daily multi vitamin traveling through your bowels and diffusing into your bloodstream These physical processes are better described by the convection diffusion equation t x L XX7 Here a is the convection speed and b is the diffusion coefficient m Convection7Diffusion 7 2234 The Convection Diffusion Equation Numerics 2 of 3 In order to guarantee that the numerical solution of the difference scheme V31 1 7 219mm bpl 7 av 1 bM1 DOV777717 also is non increasing we must have a S 1 and by S 12 when these two conditions are satisfied we have let Vf maxm l V57 1 7 2buivi bM1 7 aiv1i bM1 anvg i S S vf 1 7 2bp bp1 7 a bp1 04 Vf n1 So that lvw l iv7 is non increasing ie the peak value of the numerical solution 535 Convection7Diffusion 1434 The Convection Diffusion Equation Numerics 3 of 3 The Convection Diffusion Equation Example 1 rnnnm dxn1 1 anm25 rn2nm dxn1 1 anm25 17mm ax mm 11411125 The condition a S 1 can be re written M22 l 8 ll 11 which is a restriction on the spatial grid spacing The quantity E corresponds to the Reynolds number in fluid flow 1 1 1 or the Peclet number in heat rnsnmnxn1nanm25 17mm nxn1nanm25 r71nnmnxn1nannn25 The quantity a sometimes 2a is often called the cell Reynolds number or the cell Peclet number 39 39 X 39 f If the grid spacing h is too large then the numerical solution 1 11 1 U cannot resolve the physics and oscillations occur These oscillations are not due to instability as long as the stability 1 1 1 criterion is satleled Of course and do nOt grow they are only a Figure ForwardTime CentralSpace Convectiondiffusion with a 10 b 01 h result of inadequate resolution m 01 gt 002 k 0002541 14 lt 12 We are stable but have not resolved the physm ConvectionrDifTusion 7 2534 ConvectionrDifTusion 7 2634 The Convection Diffusion Equation Example 2 The Convection Diffusion Equation Upwind Differences 1 of 3 r n 1 m 7 nn2 11 mm 1 7 Hum 31 7 nn2 11 mm 1 7 man 31 7 nn2 11 mm In example 2 we had to push the resolution to h 002 601 points in 7111 and k 00001 10001 timesteps in 07 for a grand total of 6010601 space time grid points That is a ridiculously high price to pay for such a simple 1D problem One way to avoid the resolution restriction is to use upwind differencing of the convection term This corresponds to a switching between quot quotquotquotquot quot39quot239 quotquot39quotquotquot quot quotquotquotquot quot39quot239 quotquot39quotquotquot quot quotquotquotquotquotquot39quot239 quotquot39quotquotquot backward differencing when a gt 0 and forward differencing when a lt 0 A eg only differencing in the direction where the hyperbolic I characteristics come from j vquot1 7 vquot Vquot 7 Vquot vquot 7 vquot vquot 7 2vquot vquot7 m m m m 1 7 m1 m b m1 m m 1 k l a l h l l a l h l 112 v12 7 1 7 21w an v1quot bwrquot11 M1 2avquot1 1 1 1 Figure ForwardTime CentralSpace Convectiondiffusion with a 10 b 01 h 002 g 002 k 00001 11 14 lt 12 We are stable and have resolved the physics 5353 ConvectionrDifTusion 7 2734 ConvectionrDifTusion 7 2334 The Convection Diffusion Equation Upwind Differences 2 of 3 The restriction h 3 27b is replaced by lal 2bp lalA S 1 which is much less restrictive when b is small and a large lf w want p 14 ie k h24 then we must have h S g 1 7 g which with a 10 and b 01 as in the previous examples is 038 19 times that of the previous restriction We have however also sacrificed the spatial second order accuracy since the first order upwind difference is first order The Convection Diffusion Equation Example 3 Ynnnnnm ax ass humans Yn21u15 ax ass humans Inazn15a ax ass humans 1 ras125m ax ass humans 1quot326315 ax ass humans rauaaaaa ax ass humans 1 1 1 Figure Upwinding Convectiondiffusion with a 10 b 01 h 035 g 038 k m 0030625 11 14 lt 12 We are stable and have resolved the physics ConvectionrDifTusion 7 2934 ConvectionrDifTusion 7 3034 The Convection Diffusion Equation Example 4 The Convection Diffusion Equation Upwind Differences 3 of 3 r aamaau m am h amuaaa r ammau m am h amuaaa J J Ml ul39 M 1 1 1 Figure Upwinding Convectiondiffusion with a 10 b 01 h 040 2 038 k 004 11 14 lt 12 We are stable but have not resolved the physics 1 a2amau m am h amuaaa ConvectionrDifTusion 7 3134 The upwind scheme n1 n n 7 n Vm T V V m71 mam b n n n Vm1 T 2Vm T Vmil h2 can be rewritten in the form Vrrrl1 Vrnn Vrrrli1 V5171 7 ah Vrnn1 2V l V5171 k 3 2h 7 2 h2 We see that upwinding corresponds to changing the diffusion coefficient or adding artificial viscosity to suppress oscillations There has been much debate regarding the value of these artificial viscosity solutions clearly they may only give qualitative information about the true solution More details on solving the convection diffusion equation numerically can be found in KW MORTON Numerical Solution of Convection Diffusion Problems Chapman amp Hall London 1996 ConvectionrDifTusion 7 3234 Variable Coefficients Looking Ahead When the diffusivity b is a function of time and space eg of the common form Systems of PDEs in Higher Dimensions ut btxuxx Second Order Equations Analysis of Well Posed and Stable Problem the difference schemes must be chosen to maintain consistency For example the forward time central space scheme for this problem is given by Convergence Estimates for lVPs Well Posed and Stable lBVPs Elliptical PDEs and Difference Schemes V3 Vrnn bltnvmeZXVrnnH Vrrrli btquot7xm12Vrrrli V571 k h2 Linear lterative Methods This scheme is consistent if 1 The Method of Steepest Descent and the Conjugate Gradient btX4 S 5 Method for all values of t7X in the domain of computation m m Variable Coefficients 7 3334 Variable Coefficients 7 3434 Numerical Solutions to PDEs Lecture Notes 1 Introduction Peter Blomgren ltblomgrenpeter gmail com Department of Mathematics and Statistics ynamlcal Systems r0 Computational Science Rsearch Center San Diego State University San Diego CA 921827720 httpterminussdsuedu Spring 2010 Peter Blomgren Lecture Notes 1 7 introduction Academic Life 84 Reseach Interests 1 i i1 0 1998 2002 Research Associate Stanford University Department of Mathematics Main Focus Time Reversal and Imaging in Random Media with George Papanicolaou et al 0 1994 1998 PhD UCLA Department of Mathematics Adviser Tony F Chan PDE Based Methods for Image Processing Thesis title quotTotal Variation Methods for Restoration of Vector Valued Images quot o 1989 1994 MSc Engineering Physics Royal Institute of Technol ogy KTH Stockholm Sweden Thesis Advisers Michael Benedicks Department of Mathematics KTH and Erik Aurell Stockholm Univer sity Department of Mathematics Thesis Topic quotA Renormalization Technique for Families With Flat Maxmaquot SR Peter Blomgren Lecture Notes 1 7 introduction Outline 0 The Professor 0 Academic Life 84 Reseach Interests 9 Contact Information Office Hours a The Class Overview 0 Literature 84 Syllabus 0 Grading 0 CSU Employee Furloughs 0 Expectations and Procedures 0 Formal Prerequisites 0 Introduction 0 Some PDE Models Peter Blomgren i Lecture Notes 1 7 introduction The Professor s 7 W Academic Life t Reseach Interests ii idi Im Introduction Who is this Guy LA quot SAN Dn GD smrr m H lJMvmsiTv 9 August 2007 present Associate Professor San Diego State University Department of Mathematics and Statistics 9 August 2002 August 2007 Assistant Professor San Diego State University Department of Mathematics and Statistics Peter Blomgren i Lecture Notes 1 7 introduction The Professor Thu Intro In mm m mi I 0 Image Processing 0 Image Reconstruction 0 Noise Removal 0 Feature Extraction o Biomedical Imaging Mitochondria Heartcells 0 Time Reversal o Underwater Acoustics 0 Communications and Imaging a Wireless Communications 0 Pattern Formation a Flames in Cylindrical Burners Peter Blomgren II whim 39 7 lnM um lli m Pv l hzl39iulIIE li39ylkmht The Class 7 Overview e 552693 Finite Difference Schemes and Partial Differential Equa tions 2nd Edition Author John C Strikwerda Publisher Society for Industrial and Applied Mathematics John c sirikwcnlx Peter Blom gren The Professor I My f 39 anaemia L39xmes I Info maim Peter Bio 4am m WM 7 mmmm w 7 i 14 1 Viyn lllalglu i Chapter Title U H perbolic t rentiai Eqn 2 Analysis of Finite Difference Schemes 3 Order of Accuracy of Finite Difference Schemes 4 Stability for Multistep Methods 5 Dissipation and Dispersion 6 lic Par 39 fevers Equations 7 Systems of PDEs in Higher Dimensions 8 Second Order Equations 9 Analysis of Well Posed and Stable Problems 10 Convergence Estimates for Initial Value Problems 11 Well Posed and Stable Initial Boundary Value Problems 12 Elliptic Pattie uations and Difference Schemes 13 Linear Iterative Methods 14 Steepest Descent and Conjugate Gradient Methods l 3 Q Peter Blomgre mine We u 7 InjKtiltsdi Tire Prniessnr i i The Class 7 Overview Grading TlieEl quotw ii lrlLTI39JIlH Grading Homeworkquot 40 ProjectFinalX 60 Due almost every Friday at Noon in GMCS 587 Peter s office X Details to be determined Peter Blomgren ijji Lecture Notes 1 7 Introduction Expectations Emil Procedures Expectations and Procedures l Th F 5w V39 ilu The Class 7 Overview an Tll39 csu Employee Furloughs inin 7 CSU Employee Furloughs AY 2009 2010 Due to extraordinary budget cuts to the CSU student fees have increased 32 many sections have been cut and faculty will be required to take nine 9 unpaid furlough days each semester This is the result of a dramatic cut to the CSU by the state after years of under funding the system This class will not meet and I will not be available for office hours phone or email consultation on the following regularly scheduled days F 122 F 129 M 28 M 215 Tu 316 F 326 F 423 Th 429 Tu 54 Peter Blomgren r ji Lecture Notes 1 7 Introduction Expectations and Procedures ll 0 Most class attendance is OPTIONAL Homework and announcements will be posted on the class web page lfwhen you attend class 0 Please be on time 0 Please pay attention 0 Please turn off mobile phones 9 Please be courteous to other students and the instructor 0 Abide by university statutes and all applicable local state and federal laws Lecture Notes 1 7 Introduction Peter Blomgren l 0 Please turn in assignments on time The instructor reserves the right not to accept late assignments The instructor will make special arrangements for students with documented learning disabilities and will try to make accommodations for other unforeseen circumstances eg illness personalfamily crises etc in a way that is fair to all students enrolled in the class Please contact the instructor EARLY regarding special circumstances 0 Students are expected and encouraged to ask questions in class 0 Students are expected and encouraged to to make use of office hours If you cannot make it to the scheduled office hours contact the instructor to schedule an appointment m Peter Blomgren 3 ji Lecture Notes 1 7 Introduction 7 1225 AssJr Th The Class 7 Overview The Class The 11 III I t lntr J Expectations and Procedures In Expectations and Procedures lll Honesty Pledges l o The following Honesty Pledge must be included in all 0 Missed midterm exams Don t miss exams The instructor reserves the right to schedule make up exams make such Programs you Sme39t as Part of homework andor Prolects examsoral presentation andor base the grade solely on other g 391 your name pledge that this program is completely my own work 39nClUd39ng the nal exam work and that I did not take borrow or steal code from any a M d f I D t th f H C t Ct th other person and that I did not allow any other person to use Isse Ina exam on mIss e Ina on a e have borrow or steal portions of my code I understand that if I violate this honesty pledge I am subject to disciplinary action 355lgned pursuant to the appropriate sections of the San Diego State University Policies instructor ASAP or a grade of incomplete or F will be 0 Academic honesty submit your own work but feel free to discuss homework with other students in the class 0 Work missing the honesty Pledge may not be graded Peter Blomgren a 3 jI Lecture Notes 1 7 Introduction Peter Blomgren I39 I Lecture Notes 1 7 Introduction 39 v t1 ll 1 Resources s Expectations and Procedures Honesty Pledges ll Computer Resources You need access to a computing environment in which to write 0 Larger reports must contain the following text your code you may want to use a combination of Matlab for quick prototyping and short homework assignments and CC a your name pledge that this report Is completely my own or Fortran or something completely different work and that I did not take borrow or steal any portions from any other person Any and all references used are Class accounts for the computer labs will be available clearly cited in the text I understand that if I violate this honesty pledge I am subject to disciplinary action pursuant to the appropriate sections of the San Diego State University capable SYStem httpwwwr0hanSdsueduraCCts htmll Policies Your signature You can also use the Rohan Sun Enterprise system or another Free CC goo and Fortran 77 compilers are available for LinuxUNIX 0 Work mIssmg the honesty pledge may not be graded You may also want to consider buying the student version of Matlab http www mathworks com Peter Blomgren a 3 jI Lecture Notes 1 7 Introduction Peter Blomgren I39 I Lecture Notes 1 7 Introduction 7 1625 The Class Formal Prerequisites lriLr JI Math 693b Formal Prerequisites Graduate Bulletin Math 693b Informal Prerequisites Math 531 Math 537 and Math 693a 531 gt PDEs 0 Boundary value problems for the heat and wave equations eigenfunction expansions Sturm Liouville theory and Fourier series D39Alembert39s solution to wave equation characteristics Math 531 and Math541 or Math542 or Math543 or Math 693a Laplace39s equation maximum principles Bessel functions and Mathematical Software eg matlab 537 gt ODE u S Essential knowledge of PDEs some experience With mathematical 0 Theory of ODEs eXIstence and uniqueness dependence on H initial conditions and parameters linear systems stability and Programmmg m some language e39g39 matlab39 and near algebra39 asymptotic behavior plane autonomous systems series solutions at regular singular points Knowledge of Fourier analysis and Real analysis is not required 693a gt Advanced Numerical Analysis Numerical Optimization bUt incredibly userll 0 Numerical optimization Newton39s method for linear and nonlinear equations and unconstrained optimization Global methods nonlinear least squares integral equations 539 m Peter Blomgren t 7 3 Lecture Notes 1 7 introduction Peter Blomgren r o Lecture Notes 1 7 Introduction The Class Formal Prerequisites int Some P DE Models Introduction Some PDE Models Possibilies Finite Element Methods andi or Mimetic Finite Difference Schemes The Heat Equation Tt 7 saw T f39lt t This class will primarily focus on Finite Difference Methods The heat equation describes heat transfer in a medium H is the We will spend some time discussing Finite Element Methods thermal diffusivity and T the temperature heatmpghmovieZdric6avi andor Mimetic Finite Difference Methods in the latter part of the semester No textbook for this part handouts will be available The Wave Equation 74gt 7 V2 f x t Ci2 ff 7 The wave equation describes propagation of waves with location dependent speed wmovieZdriclavi Peter Blomgren 39 t x 3 quot Lecture Notes 1 7 Introduction Peter Blomgren r t 3 Lecture Notes 1 7 Introduction 7 Zn25 Mullah gimme WEE The Schrodinger Equation 3WX t 52 32 NET 7 w 7 VX IJX7 t i O The Schrodinger equation here in time dependent one dimensional form is the fundamental equation of physics for de scribing m mechanical behavior It is also often called the Schrodinger wave equation and is a partial differential equa tion that describes how the wavefunction IJX7 t of a physical system evolves over time new EHQWEM h is Planck39s constant divided by 27r VX a potential and the imaginary unit xil sliioviefnavi Peter Blomgren uwrm whimNil e inrrrm lrtrm 3mm elfl l39llllllilt l e Introx ucuon WEE W l b The Fokker Planck Equation 8 EGOrwrorvO 2 v VFGrrVi rorVo 2 kT V7 7GF7 F0 70 t V ETGF 7 F070 t The Fokker Planck equation describes stochastic evolution de scribing drift and diffusion of a density function G is the proba bility density F and F0 positions 7 and 70 velocities fokkerrplanckrdarkrmatleravi foklzerrplanckrlandaumpeg Peter Blomgren magnum mutant n 7 lrmxrlrrarum Amw mil il39ll bll Samaria lPtD E Me kdl ei The Korteweg de Vries Equation E 71097 with a h337 Thgp 37109 t 7 3372971 T 3 3X 3 37109 t 6t 2 3x 3 103 009 t 3X3 r The Korteweg de Vries equation governs lc 39 h is the channel height T is the surface tension g the gravitational acceleration and p the density quotPermitl nonlinear New we a The more commonly seen nondimensionalized version of the KdV equation takes the form utuxxxi6uux0 kdvstoliLonmpgkdvrs npg Peter Blomgren Exam r mum 7 Irll39f ll llHLw 1w artr cw W39azmdh Introx ucuon wx l fl DIE metals The Kuramoto Sivashinsky Equation 1 L11 V4u Vzu ElvulZ o The Kuramoto Sivashinsky equation describes the pattern forma tion of cellular flames stabilized on a circular porous plug burner mmx sign 7 innMutiltswlr uu Peter Blomgren Lemme Some PDE Models The CI Introduction Questions Comments Administrative Stuff Add Codes Available Feb 24 2nd Last day to drop classes 4th last day to add classes or change grading basis No schedule adjustments allowed after 600 pm on this date Questions Peter Blomgren Lecture Notes 1 7 Introduction Numerical Solutions to PDEs Lecture Notes 3 Hyperbolic Partial Differential Equations Convergence Consistency and Stability the CFL Condition Peter Blomgren blomgren peter gmail com Department of Mathematics and Statistics ynamical Systems roup San Diego CA 9218277720 httpterminussdsuedu Spring 2010 m Convergence Consistency and Stability 7 129 Outline 0 Recap a Convergence and Consistency 0 Introduction 0 Examples Leapfrog and LaxFriedrichs Schemes 0 The Road To Convergence 9 Stability 0 Definitions 0 The LaxRichtmyer Equivalence Theorem 0 The Courant Friedrichs Lewy Stability Condition Convergence Consistency and Stability 7 229 Last Time Classification of PDEs Hyperbolic Parabolic and Elliptic Exact solutions of the hyperbolic OneWay Wave Equation ut l ax7 tuX l bu cfx7 t constantvariable coefficient lower order term forcing Systems of hyperbolic PDEs propagation along characteristics initial and boundary values Introduction to Finite Difference Schemes gridding forward backward and central differences notation vrrr39fl 1 l avr n 7 avr n1 kh AtAX Numerical example the leapfrog scheme dependence on Convergence Consistency and Stability 7 329 Convergence and Consistency Introduction The most basic property a finite difference scheme must have in order to be useful is that Property The approximate solution to the corresponding PDE must improve as the grid spacing h7 k a 0 For now we consider inear PDEs in the form P8t3xu ft 7 X which are first order in the t derivative Also we assume that specifying the initial condition u0x u0x VX E R determines a unique solution Infinite domain no boundaries m Convergence Consistency and Stability 7 429 Convergence Definition Definition Convergent Scheme A onestep finite difference scheme approximating a PDE is a convergent scheme if for any solution to the PDE utx and solutions to the finite difference scheme vrr39n such that v converges to the initial condition u0x as m h converges to X then vrrr39 converges to utx as n k7 m h converges to tx as h7 k converge to 0 This definition is not quite complete until we clarify the convergence of vrrr39 on the discrete grid to utx which is continuously varying Convergence Consistency and Stability 7 529 Two Schemes Leapfrog and Lax Friedrichs Vr ll Vrrrliil aVr21 Vrnne1 0 2k 2h 1 1 Vrrrlll 5Vrrr39i1 V271 aVrgi1 Vrrr39i71 0 k 2h leapfrog LaxeFriedrichs Note that the Lax Friedrichs scheme is a onestep scheme whereas the leapfrog scheme is a 2 step scheme For n step schemes the definition must be changed to allow for initialization of the first n time levels Before we can apply an n step scheme we must define VP vrr39n l Vm Convergence Consistency and Stability 7 529 Grid l Initialization Suitable for 1 Step Schemes Figure Once the first level V is defined using the initial conditions v21 uoXm on the red gridepoints a onestep scheme such as LaxeFriedrichs V1711 l V17171 lav7117 k 2h can be applied to compute the values V at the remaining gridepoints one level at a time Convergence Consistency and Stability 7 729 Grid l Initialization Suitable for 2 Step Schemes Figure Once the first two levels v9 and min are defined using the initial conditions v9 uoXm and some other clever initialization scheme for v on the red grid points a twostep scheme such as the Leapfrog scheme Vn17 Vn71 Vn 7 Vni m m 3 m1 m 1 0 2k 2h can be applied to compute the values v at the remaining gridepoints one level at a time 5351 Convergence Consistency and Stability 7 329 The Lax Friedrichs Solution of the OneWay Wave Equation Lax Friednchs A0E 29 Figure Three numerical solutions of the oneway wave equation 17x Xlt1 Howl u as h k 4t 0 the solution seems to a roach the exact solution PP 535 Convergence Consistency and Stability 7 929 The Road To Convergence Directly proving that a given scheme is convergent is not easy in general However we can get there by checking two other properties consistency and stability Definition Consistent Scheme Given a PDE P8t3xu f and a finite difference scheme Pighv f we say that the finite difference scheme is consistent with the PDE if for any smooth function igtt 7 X P8t78x Pkh H 0 as k7 h gt 07 the convergence being pointWse convergence at each point 1 7 X 5351 Convergence Consistency and Stability 7 1029 Consistency Notes o For some schemes we may have to restrict how k7 h a O in order to be consistent A smooth function is a function which is differentiable at least as many times as required for the expression to make sense The difference operator PM when applied to a function of 1 7 X does not need to be restricted to grid points A forward space difference operator applied on CD at tx gives dgttx h 7 dgttx h Convergence Consistency and Stability 7 1129 Checking Consistency Lax Friedrichs 1 of 2 The Lax Friedrichs difference operator is given by gn1 fn1 1271 fn1 1271 k PLF a W 2h We use Taylor expansion around the point tmxm quotmil nm i hdgtx h2 xx i h3 m O W dgtfn1 dgt n kdgtt k2lt1gttt k3 m 0 k4 Noting that n1 1271 bin h2 xx O 74 CD 7 CD m12h e h2 m 0 hi 5351 Convergence Consistency and Stability 7 1229 Checking Consistency Lax Friedrichs 2 of 2 We can now write if 1 LF 7 7 7 pm i ta xi2k 2 k 4 dgtxx ah2 m0 74 h k2 Hence 1 h2 1 7 LF if Pdgt Pkhltigt 2kltigttt 2 k 1 2 4 74 2 bxxigah WO h Ik As long as h7 k a O in such a way that also 77 a 0 we have P 7 BEFth a 0 e the Lax Friedrichs scheme is consistent D 5351 Convergence Consistency and Stability 7 1329 Consistency 75gt Convergence 0 Consistency implies that the solution of the PDE fit is smooth is an approximate solution of the finite difference scheme FDS 0 Convergence means that a solution of the FDS approximates a solution of the PDE o It turns out that consistency is necessary but not sufficient for a FDS to be convergent We illustrate this with an example Convergence Consistency and Stability 7 1429 Example Consistency 75gt Convergence 1 of 4 We consider the oneway wave equation with constant a 1 propagation speed and apply the forward space forward time scheme 1 w 7 v Vin v 039 h A quick Taylor expansion shows that this indeed is consistent with the PDE with an error term Pd Pkh k tth xx0k2 72 We rewrite the scheme using kh n1 n k n n n n Vrn Vmigvm17 Vrn 1 Vm Vm139 Convergence Consistency and Stability 7 1529 Example Consistency 75gt Convergence 2 of 4 Let the initial condition be given by Hoof 1 if 71 x 07 0 7 O elsewhere Hence the exact solution is utx u0xi t ie a hump of height and width one traveling to the right with speed one The initial data for the FDS are given by of 1 if 71 mh 07 Vm 7 O elsewhere Convergence Consistency and Stability 7 1529 Example Consistency 75gt Convergence 3 of 4 Figure Illustration of how the exact solution propagates it is one in the band and zero outside the band The initial condition for the FDS are zeros everywhere except in the four points v8 v31 v32 and v33 where it is one 5351 Convergence Consistency and Stability 7 1729 Example Consistency 75gt Convergence 4 of 4 Figure Illustration of how the FDS solution propagates In particular we have that v i 0 n 2 0 Hence v 7Q utnxm for tnXm in the part of the band strictly in the right half plane 7 u is one there but v is zero no matter how much we refine the grid Convergence Consistency and Stability 7 1329 Stability The Missing Property If a scheme is convergent then as vrrr39 converges to utx then vrrr39 is bounded in some sense this is the essence of stability For almost all schemes there are restrictions on the way h and k can be chosen so that the particular scheme is stable A stability region is any bounded non empty region of the first quadrant of R2 that has the origin as an accumulation point Convergence Consistency and Stability 7 1929 Stability Definition Definition Stable Scheme A finite difference scheme thvr O for a first order equation is stable in a stability region if there is an integer J such that for any positive time T there is a constant CT such that 00 J 00 h E Wka Z W2 m7oo j0m00 for 0 nk T with k7 h E A Convergence Consistency and Stability 7 2029 Stability Notation The quantity 0 12 llWllh 7 Z lelZl is the L2 norm of the grid function W and is a measure of the size energy of the solution The multiplication by h is needed so that the norm is not sensitive to grid refinements the number of points increase as h a 0 With this notation the inequality in the definition can be written J 12 J 2 llvnllh CTleV lln gt llVquotllh Cle llh j0 j0 The inequality expresses a limit in terms of energy of how much the solution can grow Typically J n 71 for an n step schemem Convergence Consistency and Stability 7 2129 Checking for Stability Using the Definition Checking whether llvnllh C 20 llvjllh holds for a particular scheme directly from the definition can be a formidable task Example 151 in Strikwerda performs this test for the forward time fonNard space scheme and takes up a good page of algebraic manipulations We will return to this issue very soon with better tools in hand We note that there is a strong relation between the Stability of Finite Difference Schemes and the Well Posedness of PDEs IVPs Convergence Consistency and Stability 7 2229 Well Posedness for the NP Definition Definition Well Posed lVP The initial value problem for the first order partial differential equation Pu O is well posed if for any time T 2 0 there exists a constant CT such that any solution ut 7 X satisfies lutxl2dx CT lu0xl2dx forOgtg T All these concepts consistency well posedness stability and convergence come together in the Lax Richtmyer equivalence theorem Convergence Consistency and Stability 7 2329 The Lax Richtmyer Equivalence Theorem Theorem The Lax Richtmyer Equivalence Theorem A consistent finite difference scheme for a partial differential equation for Which the initial value problem is well posed is convergent if and only if it is stable This theorem is extremely useful Checking consistency is straight forward Taylor expansions We are going to introduce tools based on Fourier transforms which make checking stability quite easy 0 Thus if the problem is well posed we get the more difficult and desirable result Convergence by checking two relatively easy things consistency and stability 0 The relationship is one to one hence only consistent and stable schemes need to be considered Convergence Consistency and Stability 7 2429 Condition for Stability We now turn our attention to the key stability criterion for hyperbolic PDEs In last lecture we saw some numerical evidence of the leapfrog scheme applied to ut l aux 0 a 1 breaking down when gt 1 The condition lal lt 1 is necessary for stability of many explicit FDS An explicit scheme is a scheme that can be written as finite n1 n vm 7 vm n gn lmplicit schemes where the sum may contain terms with n n l 1 will be discussed soon 5351 Convergence Consistency and Stability 7 2529 The Coura nt Friedrichs Lewy Condition The following result covers all onestep schemes we have seen so far Theorem The CFL Condition For an explicit scheme for the hyperbolic equation u aux O of the form vi av21 M Mei With A kh held constant a necessary condition for stability is the Courant Friedrichs Lewy CFL condition laAl g 1 For systems of equations for Which 7 is a vector and a B and 7 are matrices we must have laiAl g 1 for all eigenvalues a of the matrix A Convergence Consistency and Stability 7 2529 The CFL Condition Proof By Picture a075 Figure Illustration of the CFL condition with A 1 held fixed The green triangle shows the region of dependence is what region influences v27 actually only the three points at the base of the triangle contribute The blue arrow corresponds to a charace teristic with speed a 075 which carries information inside the region of dependence e red arrow corresponding to a characteristic speed of a 15 carries information from outside the region of dependence 7 this information cannot be captured by the scheme m Convergence Consistency and Stability 7 2729 Another Theorem Courant Friedrichs and Lewy also proved the following theorem Theorem There are no explicit unconditionally stable consistent finite difference schemes for hyperbolic systems of partial differential equations One way of thinking about these theorems is to define the numerical speed of propagation as 1 hk and note that a necessary condition for the stability of a scheme is A71 2 lal This guarantees that the FDS can propagate information energy at least as fast as the PDE m Convergence Consistency and Stability 7 2329 Homework 1 Due 2122010 1234pm GMCS 587 Strikwerda 131 Numerical Strikwerda 141 Theoretical Strikwerda 151 Theoretical Convergence Consistency and Stability 7 2929 m r Peter Blomgren blomgren peter gmail com Department of Mathematics and Statistics ynamicai Systems roup Computationai Sciences Research Center San Diego State University San Diego CA 9218277720 harp39terminussdsuecu Spring 2010 swam Outline 0 Recap 0 Previously a The Alternating Direction Implicit Method 0 Introduction 0 CrankNicolson ADI on a 2D Square a ADI Algorithms 0 PeacemanRachford O D39Yakonov 0 Boundary Conditions for ADI Schemes 0 Implementing ADI Methods 0 PeacemanRachford 0 The MitcheIIFairweather Scheme 0 Mixed qu Derivative Terms Peter Blomgren Systems of PD The ADI Method 7 221 Recap Previously Previously We started looking at multidimensional hyperbolic and parabolic problems first via vectorvalued problems with one time and one space dimension and then to full multispace dimensional problems In terms of definitions nothing much changed the concepts of convergence consistency stability and order of accuracy are the same However some of the analysis becomes quite challenging For instance we end up needing to bound nth powers of amplification matrices llGnll g CT In order to be able to say anything useful we have to make simplifying I 5g We looked at time split schemes as a practical way to route around some of the computational challenges Stability and Boundary Conditions are a different story Peter Blon ren The ADI Method The Alternating Direction cl mpli i am A g AI Inmlome I The Alternating Direction Implicit Method The Alternating Direction Implicit ADI method is particularly useful for solving parabolic equations on rectangular domains but can be generalized to other situations Given a parabolic equation ut V o BVU 7 I711 I712 8x 7 Ht 7 axay blz by 8y LI 7 b11LIXX 2b12LIXy bQQUyy for which b11b22 gt O and bgz lt b11b22 for parabolicity Initially we will consider the case bu 0 no mixed derivative on a square domain Peter Blomgren Systems of PDEs he ADI Method 7 421 u p Direction Implicit Method H U 339 CrankrNIcolson ADI on a 20 Square Implementin Crank Nicolson on a Square If we use the Crank Nicolson schemes for 2 spatial dimensions we end up having to invert a penta diagonal matrix in each iteration in an in an in an Zn In 2n m 156 m 221 m 221 Figure LEFT The matrix which must be inverted in each CrankiNicolson iteration If we trade storage of the LUifactorization CENTER RIGHT for speed then here with 6 X 6 interior points we end up needing more than 4 times the storage For 100 X 100 interior points the requirement jumps from 49600 matrix entries tojust over 2000000 a factor of 40 The bandwidth grows linearly in n and the LUifactorization fills in the whole bandwidth In 3D the story gets even worse 7 with n X n X n interior points the bandwidth is n2 of PDE in The ADI Method n i CrankrNicolson ADI on a 2D Square The Alternating Direction Implicit Method DI em 3 Im menuquot AI The ADI Method on a Square The ADI method reduces an n dimensional problem to a sequence of n one dimensional problems We here present the idea in 2D Let A1 and A2 be two linear operators eg 82 82 Aubiu Aubiu 1 18X2 7 2 zayz For the argument to make sense we must require that we have efficient convenient ways of solving the equations Wt AW7 172 With A1 and A2 as above and a Crank Nicolson step these solutions are given by inversion of tri diagonal matrices The ADI Method The Alternating Direction Implicit Melh Dl A od hmle quot 5 son ADI on a 2D Square The ADI Method on a Square The ADI method will give us a way to solve the combined equation ut Alu l Agu7 using the available 1D solvers as building blocks Crank Nicolson applied to the combined equation gives us un1 7 Mn 1 1 k A1Lln1 Alunl E Agunll Agun O Which with some rearrangement can be written k k I 7 5A1 7 5A2 W1 1 gm 142 u 0 k3 Peter Blomgren Systems of PD The ADI Method 7 721 The Alternating Direction lmpli 39 r DI m i D Implumemi r Cranky colson ADI on a 2D Square The ADI Method on a Square Now we notice that 1i 311i 321i 31 i 32 3132 By adding and subtracting k2A1A2u on both sides of the Crank Nicolson expression we get k k k2 I 7 EA1 7 3A2 IA1A2 W1 k k k2 5A1 5A2 T141142 Lln k2 1 3 TA1A2LI 5351 Peter Blomgren Systems of PDES in D The ADI Method 7 321 The Alternating Direction lmpl r I m i I Crank colson ADI on a 2D Square mplomen The ADI Method on a Square We can factor this and use the fact that u 1 u l O k to embed the last term on the right hand side into the O k3 term k k I 7 5A1 I 7 5A2 Now if we want to advance the solution numerically we can discretize this equation and here when A1 bluxx A2 bguyy the matrices corresponding to I 7 k2 A will be tridiagonal and can be inverted quickly using the Thomas algorithm n k k n u1l A1I A2u 008 We get the discretized ADI scheme k k anl I l EALh I EALh V m k k I 7 5AM I 7 EA Peter Blomgrel Systems of PDES in The ADI Method 7 921 PeacelnawRachford Al ilhms Implementing ADI Methods ADI Algorithms Peacema n Rachford There are several approaches to solving the ADI scheme one commonly used approach is the Peaceman Rachford algorithm which also explain the origin of the name alternating direction implicit method k k IiiALh WV HEAle IPSAM v 1 II AM WW In the first half step the X direction is implicit and the y direction explicit and in the second half step the roles are reversed Is this scheme equivalent to the ADI scheme we derived It looks quite different m Peter Blomgren Systems of PD The ADI Method 7 nu21 lam ADI Algorithms Peacema n Rachford We have k I 7 5AM v 12 I I SAM v I 7 SAM WW k V 1 5A1 Hence k k I 7 5AM I 7 EA k k 7 I 5AM I 7 5AM Vn12 k k V 1 I 7 EALh I EALh WV I SAM I gAM v Note that we do not need ALhAzh ANAL1 for this to hold m Peter Blomgren Systems of PDES in The ADI Method 7 1121 The Alternating Di ion m 39 A lmpl um ADI Algorithms D39Ya konov The D39Yakonov scheme is a direct splitting of the ADI scheme we originally derived k I 7 5AM WV 19th lw ml vquot Pig42h Vn1 Vn12 Other ADI type schemes can be derived starting with other basic schemes we worked from Crank Nicolson eg the Douglas Rachford method Strikwerda pp 175 176 is derived based on backward time central space Peter Blomgren Systems of PDES in D The ADI Method 7 1221 ons for ADI Schemes The Alternating 39 ection lumlicit wk ADI A Boundary Conditions for ADI Schemes Here we consider Dirichlet boundary conditions u txy specified at the boundary in the context of the PeacemanRachford scheme k k 7 EAM WV MEAN v I 7 SAN W1 I SAM WV The correct boundary conditions for the halfstep quantity is given by 1 k 1 k n12 7 7 n 7 7 7 n1 V 2 I 24211 5 2 I 224m 5 Where did that come from Flip the second equation in the scheme add the two and solve for v 12 And it makes sense half the condition comes from the past and half from the future m Peter Blomgren Systems of PDES he ADI Method 7 1321 The mammung Direction lmpli 39 r Implementng AD Implementing ADI Methods We consider PeacemanRachford on a grid where bem ZAX7mAy 07Hi7L m OpiqM We let M kAXZ My kAyz Further we let viquot denote the fullstep quantity and W1quot denote the halfstep quantity if we are not interested in saving the results for all t kn we can overwrite these quantities We get the first halfstage b X b x i l 1 W71m l 1 MIX Wzm l 1 Wllm b 122 V171quot 1 law V1quot l 22 Vi1m forZ1HiL71andm1uiM71 Peter Blomgren Systems of PDES he ADI Method 7 1421 ion limpiicit Metl o PeacelnawRachford 20 mi Implementng AD Implementing ADI Methods Figure Active points in the first halfestep the interior points are active both for the old velayer and the welayer which is being computed Also the boundary values at the top VLM and bottom vgwo boundaries are active and so are wow 1 left and WLV right m The ADI Method PeacelnawRachford i i Implementing ADI Methods If we enumerate our grid points in the following way then we get M 71 tridiagonal systems one for each row with Lil unknowns m of PDE in The ADI Method Implementing ADI Methods We also need the missing boundary conditions for W 4 2 b n 1 n n e 5011 ml 7 Both 1 7 b W0m My 63 6am 63M 4 b2 n 1 b n b2 n M m ib 62 61m My 62M Form1HiM71mOand mMarenot needed Peter Blomgren Systems of PD The ADI Method 7 1721 The Alternating Divan 39 Implementi r Implementing ADI Methods The second half stage is given by b b i 22W VLmA 1 hwy VLm 22W Vlm1 b b 12W Wlm71 17 b1MxWzm i 12quot 7 Wlm1 for 1L71and m1M71 With the correct grid ordering we get L 7 1 tridiagonal systems of size M7 1 Boundary conditions for v are given at timelevel n i 1 Peter Blomgren Systems of PDES in D The ADI Method 7 1321 a The Alternating ion huplicit Math 7 39itluns Implementng ADI Methods 2491 it c a 0000000 00000 0000000 Figure Active points in the second halfistep left and the appropriate enumeri ation order of the gridipoints right vmaw mm 7 r eather Sch em e The M itcheII Fairweather Scheme In Strikwerda pp 180 181 there is a discussion of the MitchellFairweather scheme which is an ADI scheme which is second order in time and fourth order accurate in space 71 71 22 n12 l l 22 n 1 2ltb1uX 6gth6Xv 71472 bzuy6 h yv7 71 71 22 n1 l l 22 n12 1 2ltbzuy 6gth6yiv 71472 blux6 h6X v with Dirichlet boundary conditions for v 12 1 1 1 1 WW2 2 blMx a 1 E My 6 hz ii 3 MM 1 1 1 bl14X 7 7gt 17 7 bw 7 7gt 62 BM 6 2 y 6 y 535 Peter Blomgren Systems of PDES il D The ADI Method 7 Zn21 1 In I pl Ni 0 DI Al in n I ADI Methods Mixed Derivative Terms is Implementing ADI with Mixed uxy Derivative Terms It has been shown that no ADI scheme involving only the time levels n 1 and n can be secondorder accurate when bu y 0 ie when we have mixed derivatives A secondorder accurate modification of the PeacemanRachford scheme is given by k k 3 1 17 519116 WV 1 519226 v kblz ox oy byquot 7 Evn 1 17 szz i v 1 1 b115il WIN2 kb1260gtlt60y Vquot Vnill with Dirichlet boundary conditions for v 12 1 k 1 k Vn12 E 1 51mg n E lt17 Ebzzx n1 Peter Blomgren Systems of PD The ADI Method

### 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

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "I bought an awesome study guide, which helped me get an A in my Math 34B class this quarter!"

#### "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."

#### "Their 'Elite Notetakers' are making over $1,200/month in sales by creating high quality content that helps their classmates in a time of need."

### Refund Policy

#### STUDYSOUP CANCELLATION 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 support@studysoup.com

#### STUDYSOUP REFUND POLICY

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: support@studysoup.com

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 support@studysoup.com

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.