Page - 180 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

180 5 SolvingPartialDifferentialEquations b) TheBackwardEuler,ForwardEuler,andCrank-Nicolsonmethodscanbegiven aunified implementation. For a linearODEu0 Dau this formulation isknown as the rule: unC1 un t D .1 /aunC aunC1 : For D 0we recover the ForwardEulermethod, D 1 gives theBackward Euler scheme, and D 1=2 corresponds to theCrank-Nicolsonmethod. The approximation error in the rule is proportional to t, except for D 1=2 where it isproportional to t2. For 1=2 themethod is stable forall t. Apply the rule to theODEsystem for a one-dimensional diffusion equation. Identify the linear systemtobesolved. c) Implement the rulewithaidof theOdespypackage.Therelevantobjectname isThetaRule: solver = odespy.ThetaRule(rhs, f_is_linear=True, jac=K, theta=0.5) d) Consider thephysicalapplicationfromSect.5.1.4.Runthiscasewith the rule and D1=2 for thefollowingvaluesof t: 0.001,0.01,0.05.Reportwhatyou see. Filename:rod_ThetaRule.py. Remarks Despite the fact that theCrank-Nicolsonmethod, or the rulewith D 1=2, is theoretically more accurate than the Backward Euler and Forward Euler schemes, itmay exhibit non-physical oscillations as in the present example if the solution isvery steep. Theoscillations aredamped in time, anddecreaseswithde- creasing t. Toavoidoscillationsonemusthave t atmaximumtwicethestability limit of the Forward Euler method. This is one reasonwhy the Backward Euler method (or a 2-step backward scheme, seeExercise 5.3) are popular for diffusion equationswithabrupt initial conditions. Exercise5.6:Compute thediffusionofaGaussianpeak Solve the followingdiffusionproblem: @u @t Dˇ@ 2u @x2 ; x2 . 1;1/; t 2 .0;T (5.34) u.x;0/D 1p 2 exp x 2 2 2 ; x2 Œ 1;1 ; (5.35) @ @x u. 1;t/D0 t 2 .0;T ; (5.36) @ @x u.1;t/D0 t 2 .0;T : (5.37) Theinitialconditionis thefamousandwidelyusedGaussianfunctionwithstandard deviation (or “width”) ,which is here taken to be small, D 0:01, such that the initial condition isapeak.Thispeakwill thendiffuseandbecomelowerandwider. Computeu.x;t/untilubecomesapproximatelyconstantover thedomain. Filename:gaussian_diffusion.py.