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

220 8 SolvingOrdinaryDifferentialEquations time point tn. The second idea is that we want to progress the solution from tn to tn+1 asastraight line. We know that the line must go through the solution at tn, i.e., the point (tn,un). The differential equation tells us the slope of the line:u′(tn) = f(un,tn) = run. Thatis, thedifferentialequationgivesadirectformulaforthefurtherdirectionofthe solution curve. We can say that the differential equation expresses how the system (u)undergoeschangesat apoint. There is a general formula for a straight liney= ax+bwith slopea that goes through the point (x0,y0):y = a(x−x0)+y0. Using this formula adapted to the presentcase, andevaluating the formulafor tn+1, results in un+1 = run(tn+1 − tn)+un=un+Δtrun, which is nothing but the Forward Euler formula. You are now encouraged to do Exercise 8.2 to become more familiar with the geometric interpretation of the ForwardEulermethod. 8.2.5 ProgrammingtheFEScheme;theGeneralCase Ourpreviousprogramwas justa simplemainprogramwithout functiondefinitions, tailored to a special differential equation. When programming mathematics, it is always good to consider a (large) class of problemsand making a Python function tosolveanyproblemthatfits intotheclass.Morespecifically,wewillmakesoftware for theclassofdifferentialequationproblemsof the form u′(t)=f(u,t), u=U0, t∈[0,T], for some given functionf , and numbersU0 andT . We also take the opportunity to illustrate what is commonly called a demo function. As the name implies, the purpose of such a function is solely to demonstrate how the function works (not to be confused with a test function, which does verification by use of assert). The Python function calculating the solution must take f , U0, Δt, and T as input, find the correspondingNt, compute the solution, and return an array with u0,u1,.. .,uNt and an array with t0,t1,.. .,tNt . The Forward Euler scheme reads un+1 =un+Δtf(un,tn), n=0,.. .,Nt −1 . Thecorrespondingprogrammaynowtake the form(fileode_FE.py): import numpy as np import matplotlib.pyplot as plt def ode_FE(f, U_0, dt, T): N_t = int(round(T/dt)) u = np.zeros(N_t+1) t = np.linspace(0, N_t*dt, len(u)) u[0] = U_0