Page - 306 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

306 9 SolvingPartialDifferentialEquations b) The BackwardEuler,ForwardEuler, and Crank-Nicolsonmethodscan be given a unified implementation. For a linear ODEu′ = au this formulation is known as theθ rule: un+1 −un Δt = (1−θ)aun+θaun+1 . For θ = 0 we recover the Forward Euler method, θ = 1 gives the Backward Euler scheme, and θ = 1/2 corresponds to the Crank-Nicolson method. The approximation error in the θ rule is proportional to Δt, except for θ = 1/2 where it isproportional toΔt2. Forθ≥1/2 themethod is stable forallΔt. Applytheθ rule to theODEsystemforaone-dimensionaldiffusionequation. Identify the linear systemtobe solved. c) Implement theθ rule with aid of the Odespy package.The relevant objectname isThetaRule: solver = odespy.ThetaRule(rhs, f_is_linear=True, jac=K, theta=0.5) d) Consider the physical applicationfrom Sect.9.2.4.Run this case with theθ rule andθ = 1/2 for the followingvaluesofΔt: 0.001,0.01,0.05.Report what you see. Filename:rod_ThetaRule.py. Remarks Despite the fact that the Crank-Nicolson method, or the θ rule with θ = 1/2, is theoretically more accurate than the Backward Euler and Forward Euler schemes, it may exhibit non-physical oscillations as in the present example if the solution is very steep. The oscillations are damped in time, and decreases with decreasingΔt. To avoid oscillations one must haveΔt at maximumtwice the stability limit of the Forward Euler method. This is one reason why the Backward Euler method (or a 2-step backward scheme, see Exercise 9.3) are popular for diffusionequationswithabrupt initial conditions. Exercise9.6:Compute theDiffusion ofa GaussianPeak Solve the followingdiffusionproblem: ∂u ∂t =β∂ 2u ∂x2 , x∈ (−1,1), t∈ (0,T] (9.34) u(x,0)= 1√ 2πσ exp ( − x 2 2σ2 ) , x∈[−1,1], (9.35) ∂ ∂x u(−1,t)=0, t∈ (0,T], (9.36) ∂ ∂x u(1,t)=0, t∈ (0,T] . (9.37)