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

132 4 SolvingOrdinaryDifferentialEquations thevaluesat tn and tnC1: f.unC 1 2 ;tnC12/ 1 2 .f.un;tn/Cf.unC1;tnC1//: This results in unC1 un t D 1 2 .f.un;tn/Cf.unC1;tnC1//; which ingeneral is anonlinear algebraic equation forunC1 iff.u;t/ is not a lin- ear function ofu. To dealwith the unknown termf.unC1;tnC1/, without solving nonlinear equations, we can approximate or predict unC1 using a Forward Euler step: unC1 DunC tf.un;tn/: This reasoninggives rise to themethod u DunC tf.un;tn/; (4.57) unC1 DunC t 2 .f.un;tn/Cf.u ;tnC1//: (4.58) Theschemeapplies tobothscalar andvectorODEs. For an oscillating systemwith f D .v; !2u/ the file osc_Heun.py imple- ments thismethod.Thedemo functionin thatfile runs thesimulationfor10periods with20 timestepsperperiod. Thecorrespondingnumericalandexactsolutionsare shown in Fig. 4.21. We see that the amplitude grows, but not asmuch as for the ForwardEulermethod.However, theEuler-Cromermethod ismuchbetter! We should add that in problems where the Forward Euler method gives sat- isfactory approximations, such as growth/decay problems or the SIRmodel, the 2nd-orderRunge-KuttamethodorHeun’smethod,usuallyworksconsiderablybet- ter and produces greater accuracy for the same computational cost. It is therefore averyvaluablemethod tobeawareof, although it cannot competewith theEuler- Cromer scheme foroscillationproblems. Thederivationof theRK2/Heun scheme isalsogoodgeneral training in“numerical thinking”. 4.3.6 SoftwareforSolvingODEs There is a jungle of methods for solving ODEs, and it would be nice to have easy access to implementationsof awide rangeofmethods, especially the sophis- ticated and complicated adaptivemethods that adjust t automatically to obtain a prescribed accuracy. ThePython packageOdespy5 gives easy access to a lot of numericalmethods forODEs. The simplest possible example onusingOdespy is to solveu0 D u,u.0/ D 2, for100 timestepsuntil t D4: 5https://github.com/hplgit/odespy