Numerical Solution of Partial Differential Equations
Numerical Solution of Partial Differential Equations MA 584
Popular in Course
Popular in Mathematics (M)
This 23 page Class Notes was uploaded by Braeden Lind on Thursday October 15, 2015. The Class Notes belongs to MA 584 at North Carolina State University taught by Staff in Fall. Since its upload, it has received 11 views. For similar materials see /class/223724/ma-584-north-carolina-state-university in Mathematics (M) at North Carolina State University.
Reviews for Numerical Solution of Partial Differential Equations
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/15/15
Finite Difference Illethod 55 8 Finite difference methods for parabolic partial differential equations 81 Introduction A linear parabolic PDE has the following form m Lu where L is a linear elliptic differential operator We list some examples below I One dimensional heat equation with a source m um fct 81 The dimension is referred to the space variable even though there are two independent variables a and t I The general one dimensional second order partial differential equations the following form ac u 2bc uu Cac u lower order terms fc i If 2 an E U in the entire domain then the equation is parabolic I A general heat equation in any dimensions can be written as m Vu Vu fxt 82 where is called the heat conductivity fxt is called a source term including a sink as well I A diffusion and advection equation has the following form m V Vu WVufxt where V Vu is the diffusion term and w Va is called the advection term I A canonical form of a diffusion and reaction equation is m V Vu fxtu The non linear source term is called a reaction term I A steady state solution meaning m U of a parabolic PDE is an elliptic PDE 56 Z Li 82 Initial and boundary conditions and dynamic stability For time dependent problems7 we should have an initial condition7 usually at t 07 that is7 ux U u0x is given We also need boundary conditions as well As an example7 for one dimensional heat equation m um7 a lt a lt b we should have boundary conditions at a a and a b7 and an initial condition at t 0 There is consistency condition at 10 and 60 For example7 if a Dirichlet boundary condition is prescribed at a a and a I such that uat 91t and ub t 926 then the consistency condition is mom 97107 and u0b 92 The dynamical stability The fundamental solution for one dimensional heat equation is m an is Egg The solution is uniformly bounded However7 for the backward heat equation7 m rm if uc 0 07 then limtnoc Maxi 00 We call the solution is dynamically unstable if Maxi is not uniformly bounded In other words7 we can not find a positive constant such that 3 C There are some applications of unstable problems7 sometimes called blow up problems We will not discuss how to solve those dynamically unstable problems numerically here We will only concentrate dynamically stable problems Some commonly used finite difference methods will be discussed in this section are listed below I the forward7 backward Euler s methods I the Crank Nicolson and the 0 method a the method of line MOL if a good ODE solver can be applied a the alternating directional implicit ADI method for high dimensional problems We can use the finite difference methods for elliptic problems to take care of the spa tial discretization and the boundary conditions The crucial discussion here is the time discretization In terms of the stability of numerical methods7 we will use the Fourier transformation and the von Neumann stability analysis 83 The Euler s method Consider the heat equation m gnu fct a lt c lt b We t 2 916 no t 2 92m m 0 u0c Finite Difference Alethod 57 MOL BC C k 1 0 O O O k 0 O O 0 FW CT BWCT 10 Figure 12 A diagram of the finite difference stencils of the forward7 backward Euler method7 and the method ofquot line MOL approach We want to approximate the solution at certain time T7 0 lt T7 or all the solution between 0 lt t lt T As the first step7 we generate a grid b xizamh i01m h a m T tkzkm k01m Atzg However7 we should know that we can not use arbitrary At or n for explicit methods because of the stability concerns The second step is to approximate the derivatives with finite difference formulas We know how to discretize the spatial derivatives Let us try different finite difference formulas for the time derivative Forward Euler s method FTCT Mini tk At Mini tk At an tic Titk ulti1tk 224351 tk uci1tk h2 a u The local truncation error is 112 At TC tk uu hot 58 Z Li The discretization is order 0012 At Often we call the discretization is first order in time7 second order in space The finite difference equation then is Ur Ur UL 2U an At h where ff f an tk The solution ofquot the finite difference equations UC is an approximation 3 if 83 to the exact rib5139 tk Note that when k 07 we have the initial condition at the grid points will If we know the solution at the time level k7 the solution ofquot the finite difference equation at the next time level is Uquot 2U ltUlt3 Ut1UfAtowff i12m 1 84 Therefore we can get the solution of the finite difference equations directly om the ap proximate solution at previous time steps We do not need to solve a system of equations Such a method is called explicit time marching method Remark 81 According to our de nition the local truncation is nct At nc t nc h t 2nct nc ht At It2 However in the literature there is another de nition using nc ht 2nct nc ht h2 Tc t x t 0012 At Tct nct At nct At fct O Ath2 At The di erence is the factor of At According to this de nition the local truncation error is one order of At higher Remark 82 If fct E U and is a constant then m n and nu 613 g Wan and the local truncation error is m t Wm 1 2 0 At2 114 85 2 12 Therefore if is a constant we can choose At to get 0h4At2 0h4 0At2 order of accuracy without increase computational cornpleacity This is signi cant for the explicit method If we try the numerical method with different At and check the error against a problem that we know its exact solution7 we see the method works for some At and blows up for some other At Since the method is consistent7 it has something to do with the stability Intuitively7 we can see that to prevent the errors in om getting amplified7 we should set mm 112 lt A lt i 12 1 or t 2g 8 6 Finite Difference Illethod 59 The backward Euler s method BW CT If we use the backward finite difference formula for m at mite we get k U H Uk 2U ltU LETJ CH k121 However7 this is equivalent to the conventional expression k k k U U1quot 2 a Hill 211139 1 U131 if 8 7 At 112 39 1 39 39 The backward Euler s method is consistent The discretization error is still 0At 12 Question Can we choose At to increase the order of accuracy M1 k11 1 S However7 we can not get u with a few simple algebraic operations because all ui are coupled together We need to solve the following tri diagonal system in order to get the approximate solution at time level k 1 1 211 u of U At ff 11ng u 1 211 u U U51 mfg u 1 211 u U U51 Atff t 88 p 1 211 p US12 U514 Atfgilg p 1 2p US31 U517 At fT flill 1195 where u and ff fcitk1 Note that we can use fcitk instead of an tkH if the method is first order accurate in time Such a numerical method is called an implicit time marching method because the solution at time level k 1 are coupled together What is the advantage of the backward Euler s method It is stable for any choice of At For one dimensional problems7 the computational cost is only slightly more than the explicit Euler s method if we can use an ef cient tridiagonal solver 84 The method of line MOL If there is a good solver for ordinary differential equations systems7 we can use the method of line MOL to solve the parabolic partial differential equations Given a general parabolic equation of the form what Luct fct where L is an elliptic operator Let Lh be a finite difference operator acting on a grid an a ih We can form a semi discrete system of ordinary differential equations of the 60 Z Li following form 9U139 9t In other words we only discretize the spatial variable For the heat equation with a source I Lama 11 M u f we have L L 5 the discretize system of ODE is 9U1t i 2U1t U26 916 9t f 112 h M 01 Ui l t 212226 1mm an t 239 2 3 m 2 89 8U 716 U 720 21146 926 quotgit 5 7 fawn The initial condition is U40 u0ci0 i12 m l 810 The ODE system can be written as a vector form dy a I fyatu gt40 2 yo 811 The MOL is especially useful for non linear PDEs of the form m f mt For linear problems we typically have dy 7 A 6 dt y where A is a matrix and c is a vector Both A and c may depend on i There are many ef cient solvers for a system of ODES Most of them are based on high order Runge Kutta methods with adaptive time steps For example we can use ODE suits in Matlab and dsodej which is available through Netlz39b in Fortran It is important to know that the ODE system om the MOL is typically stiff meaning that the eigenvalues of A has very different scales For example for the heat equation the eigenvalues are between 01 and 01 12 In Matlab we can call the ODE solver using the format 1 y ode23s yfunmol 0 tfina1 yo The solution in the last row of y which can be extracted using mrnc sizey ysol ymr is the approximate solution at time t Lfinal To define the ODE system of the MOL we should create a Matlab file yfun molm whose contends contain the following Finite Difference Illethod 61 function yp yfunmolty global m h x k 1engthy ypsizek1 yp1 2y1 4 y2hh 4 ftX1 4 g1thh for i2zm2 ypi yi1 2y1 4 y2hh 4 ftXi end ypm1 ym2 2ym1 hh ftxi g2thh where glt and g2t are two matlab functions for the boundary condition at a a and a I and is the source term The initial condition can be defined as global m h x for i12m1 y0i u0xi end where u0c is a Matlab function of the initial condition 85 The CrankNicolson scheme The time step constraint At h22 for the explicit Euler s method is generally consid ered as a severe restriction If h 001 and the nal time is T 107 and 1007 we need 2107 steps The backward Euler s method does not have the time step constraint but it is only first order accurate If we want second order accuracy 0012 we need to take At 0012 Can we derive a nite difference scheme which is second order accurate both in time without compromise the stability and computational complexity The answer is the Crank Nicolson scheme The Crank Nicolson scheme is based on the following lemma which can be proved easily using the Taylor expansion Lemma 81 Let 45t be a function that has continuous rst order derivative that is 45t E 01 then A At At 2 45a 2 1 45a l w 7 Quqt hot 812 2 2 2 8 The Crank Nicolson scheme approximate the partial differential equation m fct 62 Z Li at 321 At 2 and use the averaging lemma above to approximate the spatial derivative V VuD and f t Therefore it has the following form At 2112 kl kl kl kl kl kll kl igUH 511 glUi igUi1 2m ff f1 k 39 Hf Uf 7 15 U151 Kali all We ff 813 The discretization is second order in time and second in space This can be easily proved using the following we take 1 for simplicity of the proof A quott t WutxtA2 20m4 A 3 2 ht 2 t ht W WM 0012 uc htAt 2umtAiumhtAt 2h2 uuctAt 0012 Wm uuw As mm mm 0mm QM A toll A m t fct Au 2 m t AtQ OAt2 The Crank Nicolson scheme is an implicit method In the next section7 we will prove that it is unconditional stable for the heat equation At each time step7 we need to solve a tridiagonal system of equations to get mic The computational cost is only slightly more than the explicit Euler s method7 we can take At N h and have second order accuracy It is much more ef cient than the explicit Euler s method The 0method The 0 method for the hear equation m um f t has the following form UM U k 2 k 2 k 1 k k 1 1A t1 mum 1 mum at 1 0m When 0 17 the method is the explicit Euler s method When 0 07 the method is the backward Euler s method When 0 127 the method is the Crank Nicolson scheme If the 0 3 7 then the method is unconditional stable7 otherwise7 it is conditional stable meaning there is a time step constraint The 0 method is generally first order in time and second order in space except for 0 12 Finite Difference Illethod 63 86 Stability analysis for timedependent problems Discrete Fourier trans form and vonNeumann analysis Let us first review Fourier transform FT in continuous space Let E L2 oooo 0C that is quaC lt 00 or lt 00 The Fourier transform is defined as OC A i 7 0c e iWu as as W m LC d 814 where i 1 We can is defined in the space domain while 21W is defined in the equency domain Note that if a function is defined in the domain 0 00 usually we can use the Laplace transform The inverse Fourier transform is defined as 1 0c iwcr A u as 7 e u w dw 815 V 27139 so Parseval s relation Under the Fourier transform7 we have Mb I HuH2 or f de f u2dc 816 Fourier transform is a useful tool for theoretical and numerical analysis Use the Fourier transform7 we can get rid of one derivative From the definition of Fourier transform FT7 we have 52 73 A n 8w zcu am mm 811 To show the equalities above7 we use the definition of the the inverse Fourier transform for to get 9n 7 1 0c 1W 391 V277 so 6 93 On the other hand7 since and Mad are both in L2 oo 00 we can take the partial dw derivative of the inverse Fourier transform with respect to a to get 9n 1 foo 9 1m 1 foo A l as 7 6 dw 7 zwuelm dw d E x27r 70c 05 x27r 70c du A Since the Fourier transform and its inverse transform are unique7 we have 671 iw The proof of the first equality will be left as an exercise It is easy to generalize the equality to 87 957 In other words7 we can get rid of the derivatives of one variable The Fourier transform can mm a 818 be used to study the behavior of differential equations 64 Z Li Example 1 The wave equation utauxU ooltcltoo tgtU uc0u0c If we apply the Fourier transform to the equation and the initial condition7 we get m20 or Lkaiw z Therefore we get an ODE for 22w which can be solved easily Mu t Ma U 6711quot 210W eimwquot Therefore the solution to the original wave equation is 0c neat euer mud e awtdw OC 1 foo ewmiat 2100 dw V 277 70c uc at 0 according to the definition of the inverse Fourier transform Thus we have solved the problem The solution tell us that for a pure wave equation7 the solution does not change shape but simply propagates along the characteristic line a at 0 Also note that HuH2 H lb ll wa kimwlb ll wamlb Huollz Example 2 The heat equation diffusion equation Consider m u 00 lt 15 lt 00 tgt U uc0 u0c lim u U 139 aoc We apply the Fourier transform to the equation and the initial condition to get EL 622 or m iw2 w2 The solution of the ODE then is and t 2 21m 0 air3w Therefore M2 H lb H w10e 5 2 l2 g Huollb if gt U Actually7 it can be shown that limtnoc 0 if gt 07 and limtnoc 00 if lt 0 and Hug g U Finite Difference Illethod 65 Example 3 Dispersive waves Consider 82m1 u 82m W W 101 at I where m is a non negative integer For example m at we have A m uuxa or m zw3 zw321 The solution of the ODE then is A iiw3ll uwt uw0 e Therefore Hulk ll ll2 H w10H2 HultquH2a The solution to the original PDE can be expressed as uc t an 210W e iwgl dw V 277 70c 1 foo eiwuiwm 210W dw x27r 70c It can be interpreted in the following way The Fourier component with wave number w is propagating with the velocity LL12 Waves interact with each other7 but there is no diffusion Example 4 The PDE with even order derivatives highest Consider 82m 82m 71 u 812711 835271171 m a 101 where m is a non negative integer Then we know that 2 aw2m ifm2k1 EL aiw m 2 awzm ifm2k Therefore we have A mamaaiwmhum ifm22k1 u w0 emme ifm 2k From the equality above7 we can conclude that m um and m umu are dynamically stable7 while m um and m are dynamically unstable 66 Z Li 87 Discrete Fourier transform The motivations to study a discrete Fourier transform include the stability analysis of finite difference schemes7 data analysis in equency domain7 filtering techniques etc De nition of a grid function Let quot39U721U711U01U11U21quot 39 1 be a continuous function 2135 at an i h7 the discrete Fourier transform is defined as 1 0c 235 7 Z heeimuj 819 x27r Fix Remark 83 I The de nition is d 1 J W 39 quot to the quot case that is we up proacimdte by Z and C133 by h o is continuous and periodic function of with period 27rh eiijh 27rh eiiy hgem jn eiigy h39 820 So we can focus on in the internal De nition of discrete inverse Fourier transform 1 W 14 7 WA d 821 v m ime MOS De nition of discrete Fourier and inverse Fourier transform for a nite se quence h is not involved 10101U11U21quot39 uUiWuUaUu39 We define 1 0 M Z eilgjnj Z e m 2194 822 m 9700 70 1 7r 1 A 21 m L em 5 45 823 Discrete Norm llvllh 824 it is often denoth as H21H2 as Well The Parseval s relation is also true W2 2 m WW6 20c iiithllvlf 825 h 77W Fix J h Finite Difference Illethod 67 88 De nition of stability of a nite difference scheme A finite difference scheme P s1112 U is stable in a stability region A if there is an integer J such that for any positive time T7 there is a constant CT independent of At and h such that J HVth 3 CT 2 WW 826 10 for any n that satisfies 0 g nAt g T with At h E A Remark 84 I The stability is usually independent of the source terms 2 A stable nite dijfference scheme means that the growth of the solution is at the most of a constant multiple of the sum of the norms of the solution at the rst J 1 steps 3 The stability region are all possible At and h The following theorem provide a simple way to check the stability of a finite difference scheme Theorem 81 If HVMIHh S HV39CHh is true for any h then the nite dijfference scheme is stable Proof From the condition7 we have MHz 3 V HHh S S llvlllh S HVOHu So if we take J 07 CT 17 we have the stability 89 Von Neumann stability analysis of nite difference methods The von Neumann stability analysis of a finite difference scheme can be described as the following process Discrete scheme gt discrete Fourier transform gt growth factor 95 gt stability g 1 And the simplification of the von Neumann analysis Example 1 The forward Euler method FW CT for the heat equation m n is U 2m U At U 2 U H 1 1 h L1 a u I Lhz 68 Z Li From the discrete Fourier transform7 we know 1 7rh A U 7 WW d 72 We a a k 1 7 1 1 Wk 12 Vic U 7 1 9 U d 7 19 61 U d 39 H1 72 i he 5 27 i he 5 Similarly 1 Wh A U 7 1m quotghU c d J 1 m i he 3 Plugging these relations into the nite difference scheme7 we obtain 1 7R A Hf 7 61W 1 037181 2 6191 Mad 27f 7mm On the other hand7 according to the definition7 we have k1 1 WI 39g h A k1 U 7 91 7 U d 1 m M a a The discrete Fourier transform is unique which implies 01115 2 1 WW1 2 M 015 2 mum where 95 2 1 1107 quot 2 elfquot is called the growth factor If 3 17 then WWI 3 Wk and thus H k llh S H th The nite difference scheme then is stable Let us examine now 95 1 u cos h isin h 2 cos h isin h 1 2 cos h 1 1 4p sin2 h2 31 But we need that 3 17 or 1 3 95 3 1 Note that 1s1 4p 31 411 sir12 h2 2 96 31 So if we take 1 3 1 411 then we can guarantee that 3 17 so the stability Therefore a suf cient condition for the stability of the forward Euler s method is h2 u while we can not claim what will happen if the condition is violated7 it is a reasonable good upper bound for the stability 131 411 or 4 3 2 or At 3 827 Finite Difference Illethod 69 Simpli cation of von Neumann stability analysis for one step time march ing method Assume that we have a one step time marching method UM1 fUkUk1 We have the following theorem that tells the stability of a finite difference method Theorem 82 Let 0 h A onestep nite di erence scheme with constant coef cients is stable if and only if there is a constant K independent of At and h and some positive grid spacing Atg and ho such that g0Ath 3 1 KAt 828 for all 0 0 lt k 3 kg and 0 lt h 3 ho Ifg0Ath is independent ofh and At the stability condition 828 can be replaced with 19W S 1 829 This theorem shows that to determine the stability of a nite di erence scheme we need con sider only the ampli cation factor 9015 90 This observation is due to non Neumann and because of that this analysis is usually called yon Neumann analysis We present some examples below We can follow the following steps for the von Neumann analysis 1 Set U ell Plugging the expression into the finite difference scheme 2 Express U as U g ei7h 3 Solve and determine whether or when g 1 for the stability Example 2 The stability of backward Euler s method for the heat equation m n is k1 k1 kz1 UM 7 Uk LN 1 r 1 h2 h2 39 Follow the procedure mentioned above7 we have germ euhg u pariah germ Gianna 95 Grok 1 H 949 2 alibi 96 70 Z Li We solve 95 to get 1 1103 h 2 elf 1 7 1 1 u2cosh 2 1 4p 51112015 2 96 lt1 for any h and At gt 0 Also it is obvious that 1 lt 0 3 95 Therefore7 3 17 and the backward Euler s method is unconditionally stable Example 3 The Leapfrog scheme two stage method for the heat equation m an is 1 k1 k k k U U1 UPI 2U Um39 2At 12 830 This method is unconditionally unstable To show this7 we follow the stability analysis Note that we need to use Ufil eijhgQQ We can get lt e 7hg 917 emh MWz 2 elm 617 eijhgilp si112h 2 We get a quadratic equation for 95 962 4H sums2 96 1 I 0 831 The solution are 2psi112h 2 t V4112 sin4h 2 1 One of roots is 95 2 2u 5113015 4H251n4h 2 1 whose magnitude 95 2 1 While we do not have any conclusion about the stability of the method7 the method is known to be a unstable method in the literature 810 Finite difference methods and analysis for 2D parabolic equations The general form of equations is m alul 12 uyy mt fc y t Finite Difference Alethod 71 with boundary conditions and an initial condition It can be written as m Lu f Where L is the spatial differential operator Therefore the method of line MOL can be used if one can find a good ODE solver for stiff ODE systems Note that the system is large 00712 For simplicity we discuss the heat equation u V Vu f t We assume is a constant The simplest method is the forward Euler s method U811 Hg H Harlin Ulim U514 Ut fw 4115 milky where u At 112 The method is first order in time and second order in space and it is conditionally stable The stability condition is h2 At 3 E 832 Note that factor of 4 instead of 2 To show this using von Neumann analysis with f 0 we should set Z 6139 lhm ljhy 2 af x where g 1515217 x 111 11317 US 2 game 834 Substituting these expressions into the finite difference scheme we can get g l 2 1 4p sin2 1h2 si112 2h2 g 1 where we assume hi by h for simplicity Therefore to get the stability we enforce 11 4p2 31 4 sin2 1h2 siI12 2h2 The inequality om the very left is 4At 2 lt 2 112 Alti h2 or t4 Backward Euler s method BWCT The scheme can be written as k1 k M1 M1 M1 k1 k1 U17 U17 UHJ UiHJ Uty el UMH 4111 At 112 The method is first order in time and second order in space and it is unconditionally stable Lyn 835 The coef cient matrix for the unknown Us is a block tridiagonal and it is strictly row diagonally dominant if we use the natural row ordering 72 Z Li CrankNicolson scheme CN The scheme can be written as k1 k k1 k1 k1 k1 k1 U17 Hit A UHJ Hem Uiaiil UMH 4111 2 h k1 At f 836 l 1 112271 U314 U514 Ui f l 4112 fig h2 19 39 Both the local truncation error and global error are 0 At2 112 The scheme is uncondi tionally stable for linear problems However7 we need to solve a system of equations whose coef cient matrix is a strictly row diagonally dominant and block tridiagonal matrix if we use the natural row ordering 811 The alternating directional Implicit ADI method The ADI is one special case of time splitting or fractional step methods The idea is to use implicit discretization in one direction while using explicit in another direction For the heat equation m nu uw f yt7 the ADI method is kl k1 Mi Mi U17 2 U U142 2 2 Ui12i U514 2111quot U57 flt At2 112 112 39 17 a 8 3H l k M kl kl kl k 1 k 1 k 1 Uii U17 2 U142 2 2 Ui12i U571 2U Uni flt At2 12 m 19 39 The method is second order in time and space is uc y t E C4 in space It is unconditionally stable for linear problems We can use symbolic expressions to discuss the method The method can be written as kl At Ml At At kl Hi 2 2 US 753 2 355M271 317 2 8 38 MI 7 1 At2 Ml At2 Ml At Ml 39 U17 U15 2 7511 2 Ydyinj fij 2 In the matrix vector form7 if we move unknowns to the left hand side7 then we get A A A I 133 Uk Li I 5135 Uk gFW 839 I D UM1 I 03 U Fk The following derivation is to get a simple form for convenience of analysis Solve for Uk Li to get L At 1 At At 1 At L U m I 71 1 71 U C I 71 F C2 Finite Difference Illethod Plugging the expression into the second equation in the ADI method to get At2 W At2 At2 1 At2 k 1 D9U 1D1 1 D1 1Dy U At At 1 At 1 At t I 7133 I 7133 7F W 713 We can go further to get At At At At 1 703 I 713 UM I 7133 I 713 Uk 2 I 133 Fmi F Note that in the derivation7 we have used the fact that At 2 At 2 At 2 At 2 I 7 I 7 I 7Dy I 7 and other communitive operations Implementation We take advantages of the tridiagonal solver M1 k At 2 M1 At 2 2 At Ml U17 2 U17 75mm 2 75WU15 17 2 k l k k l For a xed 17 we get a tridiagonal system of equations for U1 27 U1 27 7 Unify assuming a Dirichlet boundary condition at as a and as b The system of equations in the matrix vector form is Uk 1 2H H 19 k 1 U2 2 l p 12p p Ugly p 12p p p 12p p U i u 12u Ug where w3y gHJmmw H rrgU aUEH U57 fj H U574 219514 Hiya I 2 U4 f u veil 2 am At k1 U371 7fm32 H 152444 21173444 UglilfH At kl 1 k k k r Umilg 2 fij HltUm71471 21171271471 Ugim Hubcwu 311 7 5AL ML 7 kl where u W and f1 f amt 2 For each 7 we need to solve a symmetric tndlagonal system of equations Pseudo code in Matlab for j 2m 7 Look for fixed yj A sparsem1m 1 bzerosm1l for i22m bil ulijl 2ulij ulijlhl ft2xiyj 2u1ijdt if 1 bil bil uexactltt2xi1yjh1 Aili lhl else if im bi1 bi1 uexactt2xi1yjh1 Ai1i2 1h1 else Ai1i 1h1 Ai1i2 1h1 end end Ai1i1 2dt 2h1 end ut Ab Z Solve the diagonal matrix l loop in y direction Finite Difference Niethod for i A for j2n 22m sparsem1 m1 bzeros m1 1 bj1 u2i1j 2u2ij u2i1jh1 ft2xiyj 2u2ijdt bj1 uexactltt1xi yj1h1 ifj bj1 1h1 Aj1j else if jn bj1 bj1 uexactltt1xiyj1h1 1h1 Aj1j2 else Aj1j 1h1 system Aj1j2 1h1 70 Solve the end end Aj1j1 2dt 4 2111 end ut Ab 8111 Truncation error analysis If we add the two equations in 837 together7 we get UK H k l l y 2 2 k 1 k M 1 A t 2 25mm 2 5w U1 U17 2f 2 840 If we subtract the first equation onx the second equation in 8377 we get k1 k 2 k1 k 2 U17 U17 Aww Uij U17 841 l 4113 Plugging this into 8407 we can get 11 U U U E kl H J 2 5 55 2 U it 2 At 2 2 511599 Now we can see clearly that the discretization is second order accurate both in space and 1 4 time7 that is 0At2 12 842 8112 The stability analysis We take f U and set 6 139 1hu 2h2j1 Ul H M61162 eii lhll 2hj39 76 Z Li Using the operator form7 we have At At At At 1 7531 1 7559 Ugfl 1 7531 1 7559 U 1 1 559 915142 eiiglw wm At At 1 1 2 1 2 55y871 1h112h2j39 After some manipulations7 we can get 1 usiI12 1h2 1 usin2fgh2 1 usiI12 1h2 1 psit1252h2 where u and we have set It hz h for simplicity So we have g 1 2 g 1 no 01 951u52 matter what At and h are Therefore the ADI method is unconditionally stable for linear heat equations 812 An implicitexplicit mixed method for diffusion and and advection equations Consider the equation m W Vu V Vu fcyt It is not so easy to get second order implicit scheme that the coef cient matrix is diagonally dominant or symmetric positivenegative definite due to the advection term One solution is to use implicit scheme for the diffusion term and an explicit scheme for the advection term The scheme has the following form om time level t 6 to tk w W 39 thk Vh 39 W Uk Vb 39 VWWLI 4 1quot 83943 where w WM 2 g ltwvtu1k w vtu 844 The time step constraint can be chosen as h At 3 845 At each time step7 we need to solve a Helmholtz equation V Vuk1 2 ng 9 2w Vuk V Vuk 2f 846 We need 111 to get this method started We can use the explicit Euler s method FW CT to approximate 111 It should have no effect on the stability and global error 0At2 12 Finite Difference Illethod 1 1 813 Solving elliptic PDES using the numerical methods for parabolic PDES We know that the steady state solution of a parabolic PDE is the solution of the corre sponding elliptic PDE For example7 the steady state solution of the parabolic PDE m V Vu wu fxt is the solution to the elliptic PDE V Vu wufx 0 if the limit fx 2 fxu t exists The initial condition is irrevalent to the steady state solution But the boundary condition is There are some advantages of this approach especially for some non linear problems in which the solution is not unique Using this approach7 we can control the variation in the intermediate solutions The linear system of equations are more diagonally dominant Since we only care about the steady state solution7 we prefer to use implicit methods with large time steps The accuracy in time in unimportant
Are you sure you want to buy this material for
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'