Seite - 256 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Band Second Edition

256 8 SolvingOrdinaryDifferentialEquations The integralcanbeapproximatedby the famousSimpson’s rule6: tn+1∫ tn f(u(t),t)dt ≈ Δt 6 ( fn+4fn+12 +fn+1 ) . The problem with this formula is that we do not know fn+12 = f(un+12,tn+12) andfn+1 = (un+1,tn+1) as onlyun is available and onlyfn can then readily be computed. To proceed, the idea is to use variousapproximationsforfn+12 andfn+1 based onusingwell-knownschemesfor theODE in the intervals [tn,tn+12]and [tn,tn+1]. Letus split the integral into four terms: tn+1∫ tn f(u(t),t)dt≈ Δt 6 ( fn+2fˆn+12 +2f˜n+12 + f¯n+1 ) , where fˆn+12 , f˜n+12 ,and f¯n+1 areapproximationstofn+12 andfn+1 thatcanutilize already computed quantities. For fˆn+12 we can simply apply an approximation to un+12 basedona ForwardEuler stepof size 12Δt: fˆn+ 1 2 =f(un+ 1 2 Δtfn,tn+12) (8.63) This formula provides a prediction offn+12 , so we can for f˜n+12 try a Backward Eulermethod toapproximateun+12: f˜n+ 1 2 =f(un+ 1 2 Δtfˆn+ 1 2,t n+12). (8.64) With f˜n+12 as an approximation to fn+12 , we can for the final term f¯n+1 use a midpoint method (or central difference, also called a Crank-Nicolson method) to approximateun+1: f¯n+1 =f(un+Δtfˆn+12,tn+1). (8.65) WehavenowusedtheForwardandBackwardEulermethodsaswellas thecentered difference approximation in the context of Simpson’s rule. The hope is that the combination of these methods yields an overall time-stepping scheme from tn to tn+1 that is much more accurate than the individual steps which have errors proportionaltoΔt andΔt2.This is indeedtrue: thenumericalerrorgoes in fact like CΔt4 fora constantC, whichmeans that the errorapproacheszeroveryquicklyas we reduce the time step size, compared to the ForwardEuler method (error∼Δt), 6 http://en.wikipedia.org/wiki/Simpson’s_rule.