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

8.4 Oscillating1DSystems:ASecondOrderODE 249 8.4.6 SoftwareforSolvingODEs There is a jungle of methods for solving ODEs, and it would be nice to have easy access to implementations of a wide range of methods, espe- cially the sophisticated and complicated adaptive methods that adjust Δt automatically to obtain a prescribed accuracy. The Python package Odespy (https://github.com/thomasantony/odespy/tree/py36/odespy) gives easy access to a lot ofnumericalmethodsforODEs. Odespy: Example with Exponential Growth The simplest possible example on usingOdespy is to solveu′ =u,u(0)=2, for100 timestepsuntil t=4: import odespy import numpy as np import matplotlib.pyplot as plt def f(u, t): return u method = odespy.Heun # or, e.g., odespy.ForwardEuler solver = method(f) solver.set_initial_condition(2) time_points = np.linspace(0, 4, 101) u, t = solver.solve(time_points) plt.plot(t, u) plt.show() In other words, you define your right-hand side function f(u, t), initialize an Odespysolverobject, set the initialcondition,computeacollectionof timepoints where you want the solution, and ask for the solution. If you run the code, you get theexpectedplotof theexponential function(notshown). A nice feature of Odespy is that problem parameters can be arguments to the user’sf(u, t) function.Forexample, if ourODE problemisu′ =−au+b, with twoproblemparametersa andb,wemaywriteourf functionas def f(u, t, a, b): return -a*u + b Theextra,problem-dependentargumentsaandbcanbe transferredto this function if we collect their values in a list or tuple when creating the Odespy solver and use thef_argsargument: a = 2 b = 1 solver = method(f, f_args=[a, b]) This is a good feature because problem parameters must otherwise be global variables—now they can be arguments in our right-hand side function in a natural way. Exercise 8.21 asks you to make a complete implementation of this problem andplot thesolution.