8.4 Oscillating 1D Systems: A Second Order ODE
249
8.4.6 Software for Solving ODEs
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, especially 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 of numerical methods for ODEs.
Odespy: Example with Exponential Growth The simplest possible example on
using Odespy is to solve u = u, u(0) = 2, for 100 time steps until 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
Odespy solver object, set the initial condition, compute a collection of time points
where you want the solution, and ask for the solution. If you run the code, you get
the expected plot of the exponential function (not shown).
A nice feature of Odespy is that problem parameters can be arguments to the
user’s f(u, t) function. For example, if our ODE problem is u = −au + b, with
two problem parameters a and b, we may write our f function as
def f(u, t, a, b):
return -a*u + b
The extra, problem-dependent arguments a and b can be transferred to this function
if we collect their values in a list or tuple when creating the Odespy solver and use
the f_args argument:
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
and plot the solution.
249
8.4.6 Software for Solving ODEs
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, especially 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 of numerical methods for ODEs.
Odespy: Example with Exponential Growth The simplest possible example on
using Odespy is to solve u = u, u(0) = 2, for 100 time steps until 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
Odespy solver object, set the initial condition, compute a collection of time points
where you want the solution, and ask for the solution. If you run the code, you get
the expected plot of the exponential function (not shown).
A nice feature of Odespy is that problem parameters can be arguments to the
user’s f(u, t) function. For example, if our ODE problem is u = −au + b, with
two problem parameters a and b, we may write our f function as
def f(u, t, a, b):
return -a*u + b
The extra, problem-dependent arguments a and b can be transferred to this function
if we collect their values in a list or tuple when creating the Odespy solver and use
the f_args argument:
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
and plot the solution.