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

140 4 SolvingOrdinaryDifferentialEquations Westartwith integrating thegeneralODEu0 Df.u;t/overa timestep, from tn to tnC1, u.tnC1/ u.tn/D tnC1Z tn f.u.t/;t/dt : Thegoalof thecomputation isu.tnC1/ (unC1),whileu.tn/ (un) is themost recently knownvalueofu. Thechallengewith the integral is that the integrand involves the unknownubetween tn and tnC1. The integral canbeapproximatedby the famousSimpson’s rule6: tnC1Z tn f.u.t/;t/dt t 6 fnC4fnC12 CfnC1 : The problemwith this formula is that we do not know fnC 1 2 D f.unC12 ;tnC12/ andfnC1 D .unC1;tnC1/ as onlyun is available and onlyfn can then readily be computed. Toproceed, the idea is tousevariousapproximations forfnC 1 2 andfnC1based onusingwell-knownschemes for theODE in the intervals Œtn;tnC12 and Œtn;tnC1 . Letus split the integral into four terms: tnC1Z tn f.u.t/;t/dt t 6 fnC2 OfnC12 C2 QfnC12 C NfnC1 ; where OfnC12 , QfnC12 , and NfnC1 are approximations tofnC12 andfnC1 that canuti- lizealreadycomputedquantities. For OfnC12 wecansimplyapplyanapproximation tounC 1 2 basedonaForwardEuler stepof size 1 2 t: OfnC12 Df unC 1 2 tfn;tnC12 (4.63) This formula provides a prediction offnC 1 2 , sowe can for QfnC12 try aBackward Eulermethod toapproximateunC 1 2 : QfnC12 Df unC 1 2 t OfnC12 ;tnC12 : (4.64) With QfnC12 as an approximation to fnC12 , we can for the final term NfnC1 use amidpointmethod (or central difference, also called aCrank-Nicolsonmethod) to approximateunC1: NfnC1 Df.unC t OfnC12 ;tnC1/: (4.65) We have nowused the Forward andBackwardEulermethods aswell as the cen- tered difference approximation in the context of Simpson’s rule. The hope is that 6http://en.wikipedia.org/wiki/Simpson’s_rule