This notebook contains an excerpt from the Python Programming and Numerical Methods - A Guide for Engineers and Scientists, the content is also available at Berkeley Python Numerical Methods.
The copyright of the book belongs to Elsevier. We also have this interactive book online for a better learning experience. The code is released under the MIT license. If you find this content useful, please consider supporting the work on Elsevier or Amazon!
< 23.4 Numerical Error and Instability | Contents | 23.6 Summary and Problems >
Python ODE Solvers (BVP)¶
In scipy, there are also a basic solver for solving the boundary value problems, that is the scipy.integrate.solve_bvp function. The function solves a first order system of ODEs subject to two-point boundary conditions. The function construction are shown below:
CONSTRUCTION:
Let \(F\) be a function object to the function that computes
\(t\) is a one-dimensional independent variable (time), \(S(t)\) is an n-dimensional vector-valued function (state), and the \(F(t, S(t))\) defines the differential equations. \(S0\) be an initial value for \(S\). The function \(F\) must have the form \(dS = F(t, S)\), although the name does not have to be \(F\). The goal is to find the \(S(t)\) approximately satisfying the differential equations, given the initial value \(S(t0)=S0\).
The way we use the solver to solve the differential equation is: $\(solve\_ivp(fun, t\_span, s0, method = 'RK45', t\_eval=None)\)$
where \(fun\) takes in the function in the right-hand side of the system. \(t\_span\) is the interval of integration \((t0, tf)\), where \(t0\) is the start and \(tf\) is the end of the interval. \(s0\) is the initial state. There are a couple of methods that we can choose, the default is ‘RK45’, which is the explicit Runge-Kutta method of order 5(4). There are other methods you can use as well, see the end of this section for more information. \(t\_eval\) takes in the times at which to store the computed solution, and must be sorted and lie within \(t\_span\).
Let’s try one example below.
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_bvp
plt.style.use('seaborn-poster')
%matplotlib inline
F = lambda t, s: np.cos(t)
t_eval = np.arange(0, np.pi, 0.1)
sol = solve_ivp(F, [0, np.pi], [0], t_eval=t_eval)
plt.figure(figsize = (12, 4))
plt.subplot(121)
plt.plot(sol.t, sol.y[0])
plt.xlabel('t')
plt.ylabel('S(t)')
plt.subplot(122)
plt.plot(sol.t, sol.y[0] - np.sin(sol.t))
plt.xlabel('t')
plt.ylabel('S(t) - sin(t)')
plt.tight_layout()
plt.show()
def fun(t, y):
print(y[1].shape)
return np.vstack((y[1], -9.8))
f = lambda t, s: \
np.dot(np.array([[0,1],[0,-9.8/s[1]]]),s)
def bc(ya, yb):
return np.array([ya[0], yb[0]])
t = np.linspace(0, 5, 11)
y_a = np.zeros((2, t.size))
y_a[-1]= 50
res_a = solve_bvp(f, bc, t, y_a)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-23-2946ad05c3b2> in <module>
4 y_a[-1]= 50
5
----> 6 res_a = solve_bvp(f, bc, t, y_a)
~/miniconda3/lib/python3.6/site-packages/scipy/integrate/_bvp.py in solve_bvp(fun, bc, x, y, p, S, fun_jac, bc_jac, tol, max_nodes, verbose)
1053 fun, bc, fun_jac, bc_jac, k, a, S, D, dtype)
1054
-> 1055 f = fun_wrapped(x, y, p)
1056 if f.shape != y.shape:
1057 raise ValueError("`fun` return is expected to have shape {}, "
~/miniconda3/lib/python3.6/site-packages/scipy/integrate/_bvp.py in fun_p(x, y, _)
649 if k == 0:
650 def fun_p(x, y, _):
--> 651 return np.asarray(fun(x, y), dtype)
652
653 def bc_wrapped(ya, yb, _):
~/miniconda3/lib/python3.6/site-packages/numpy/core/numeric.py in asarray(a, dtype, order)
536
537 """
--> 538 return array(a, dtype, copy=False, order=order)
539
540
ValueError: setting an array element with a sequence.
def fun(x, y):
return np.vstack((y[1], -np.exp(y[0])))
def bc(ya, yb):
return np.array([ya[0], yb[0]])
x = np.linspace(0, 1, 5)
y_a = np.zeros((2, x.size))
y_b = np.zeros((2, x.size))
#y_b[0] = 3
res_a = solve_bvp(fun, bc, x, y_a)
/Users/qingkaikong/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: divide by zero encountered in true_divide
/Users/qingkaikong/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in multiply
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-27-845fe81f85d9> in <module>
----> 1 res_a = solve_bvp(F, bc, x, y_a)
~/miniconda3/lib/python3.6/site-packages/scipy/integrate/_bvp.py in solve_bvp(fun, bc, x, y, p, S, fun_jac, bc_jac, tol, max_nodes, verbose)
1053 fun, bc, fun_jac, bc_jac, k, a, S, D, dtype)
1054
-> 1055 f = fun_wrapped(x, y, p)
1056 if f.shape != y.shape:
1057 raise ValueError("`fun` return is expected to have shape {}, "
~/miniconda3/lib/python3.6/site-packages/scipy/integrate/_bvp.py in fun_p(x, y, _)
649 if k == 0:
650 def fun_p(x, y, _):
--> 651 return np.asarray(fun(x, y), dtype)
652
653 def bc_wrapped(ya, yb, _):
~/miniconda3/lib/python3.6/site-packages/numpy/core/numeric.py in asarray(a, dtype, order)
536
537 """
--> 538 return array(a, dtype, copy=False, order=order)
539
540
ValueError: setting an array element with a sequence.
res_b = solve_bvp(fun, bc, x, y_a)
res_a.y
message: 'The algorithm converged to the desired accuracy.'
niter: 1
p: None
rms_residuals: array([9.86500717e-05, 1.05360602e-04, 1.05360602e-04, 9.86500717e-05])
sol: <scipy.interpolate.interpolate.PPoly object at 0x1832058938>
status: 0
success: True
x: array([0. , 0.25, 0.5 , 0.75, 1. ])
y: array([[ 0.00000000e+00, 1.04784145e-01, 1.40534773e-01,
1.04784145e-01, 0.00000000e+00],
[ 5.49349275e-01, 2.84320977e-01, -1.28752555e-18,
-2.84320977e-01, -5.49349275e-01]])
yp: array([[ 5.49349275e-01, 2.84320977e-01, -1.28752555e-18,
-2.84320977e-01, -5.49349275e-01],
[-1.00000000e+00, -1.11047088e+00, -1.15088910e+00,
-1.11047088e+00, -1.00000000e+00]])
The above figure shows the corresponding numerical results. As in the previous example, the difference between the result of solve_ivp and the evaluation of the analytical solution by Python is very small in comparison to the value of the function.
EXAMPLE:
Let the state of a system be defined by \(S(t) = \left[\begin{array}{c} x(t) \\y(t) \end{array}\right]\), and let the evolution of the system be defined by the ODE
Use solve_ivp to solve this ODE for the time interval \([0, 10]\) with an initial value of \(S_0 = \left[\begin{array}{c} 1 \\1 \end{array}\right]\). Plot the solution in (\(x(t), y(t)\)).
F = lambda t, s: np.dot(np.array([[0, t**2], [-t, 0]]), s)
t_eval = np.arange(0, 10.01, 0.01)
sol = solve_ivp(F, [0, 10], [1, 1], t_eval=t_eval)
plt.figure(figsize = (12, 8))
plt.plot(sol.y.T[:, 0], sol.y.T[:, 1])
plt.xlabel('x')
plt.ylabel('y')
plt.show()
In practice, some ODEs have bad behavior known as stiffness. Loosely speaking, stiffness refers to systems that can have very sharp changes in derivative. An example of a stiff system is a bouncing ball, which suddenly changes directions when it hits the ground. Depending on the properties of the ODE you are solving and the desired level of accuracy, you might need to use different methods for solve_ivp. There are many methods that you can choose for the method argument in solve_ivp, take a look of the documentation to understand it more. As suggested by the documentation, you should use the ‘RK45’ or ‘RK23’ method for non-stiff problems and ‘Radau’ or ‘BDF’ for stiff problems. If not sure, first try to run ‘RK45’. If needs unusually many iterations, diverges, or fails, your problem is likely to be stiff and you should use ‘Radau’ or ‘BDF’. ‘LSODA’ can also be a good universal choice, but it might be somewhat less convenient to work with as it wraps old Fortran code.
< 23.4 Numerical Error and Instability | Contents | 23.6 Summary and Problems >