Python:Ordinary Differential Equations

IMPORTANT

It seems there is some weirdness with the RK45 (default) solver in solve_ivp with the default relative tolerance. It produces underdamped regions for the solution where it should not. Using the

rtol = 1e-5

kwarg seems to fix this so, for the moment, you should append that kwarg to the parameter list for all solutions. More information when I know more!

Introduction

This page, based very much on MATLAB:Ordinary Differential Equations is aimed at introducing techniques for solving initial-value problems involving ordinary differential equations using Python. Specifically, it will look at systems of the form:

\( \begin{align} \frac{dy}{dt}&=f(t, y, c) \end{align} \)

where \(y\) represents an array of dependent variables, \(t\) represents the independent variable, and \(c\) represents an array of constants. Note that although the equation above is a first-order differential equation, many higher-order equations can be re-written to satisfy the form above.

In addition, the examples on this page will assume that the initial values of the variables in \(y\) are known - this is what makes these kinds of problems initial value problems (as opposed to boundary value problems).

Solving initial value problems in Python may be done in two parts. The first will be a function that accepts the independent variable, the dependent variables, and any necessary constant parameters and returns the values for the first derivatives of each of the dependent variables. In other words, you will need to write a function that takes \(t\), \(y\), and possibly \(c\) and returns \(f(t, y, c)\). Note that the output needs to be returned as an array or a list.

The second part will use this function in concert with SciPy's ODE solver to calculate solutions over a specified time range assuming given initial conditions. Specifically it will use scipy.integrate.solve_ivp. Take a look at the help information and examples on that page before continuing here.

Simple Example

Assume you want to numerically solve:

\( \begin{align} \frac{dy(t)}{dt}=t-y(t) \end{align} \)

for some times between 0 and 15. Assuming the initial value for \(y(0)\) is 2, you can solve this with:

import numpy as np
from scipy.integrate import solve_ivp
sol = solve_ivp(lambda t, y: t-y, [0, 15], [2], rtol = 1e-5)

After this runs, sol will be an object containing 10 different items. Of these, sol.t will be the times at which the solver found values and sol.y will be a 2-D array. Each row of sol.y will be the solution to one of the dependent variables -- since this problem has a single differential equation with a single initial condition, there will only be one row. To plot the solution, use:

import matplotlib.pyplot as plt
plt.plot(sol.t, sol.y[0], 'k--s')

The [0] at the end of sol.y is critical as it pulls out the first row of values.

Note in this case that the solutions are not at evenly spaced times - specifically, they are at:

In [1]: sol.t.reshape(-1,1)
Out[1]: 
array([[ 0.        ],
       [ 0.0932467 ],
       [ 0.9624543 ],
       [ 1.8316619 ],
       [ 2.71989099],
       [ 3.84560123],
       [ 5.27819586],
       [ 7.1756939 ],
       [ 9.80447762],
       [13.62674388],
       [15.        ]])

The only guarantee with the solve_ivp program is that you will have points at the start and end of the range; if you want solutions at specific intermediate times, you can specify that with the t_eval kwarg:

sol = solve_ivp(lambda t, y: t-y, [0, 15], [2, 3], 
                t_eval=np.linspace(0, 15, 16), rtol = 1e-5)

will give solutions for times:

sol.t.reshape(-1,1)
Out[42]: 
array([[ 0.],
       [ 1.],
       [ 2.],
       [ 3.],
       [ 4.],
       [ 5.],
       [ 6.],
       [ 7.],
       [ 8.],
       [ 9.],
       [10.],
       [11.],
       [12.],
       [13.],
       [14.],
       [15.]])

The solver may have solved for several more points than this to get where it was going; once the t_eval kwarg is in place, it will only report at the times specified.

Another Example

Imagine you want to look at the value of some output \(y(t)\) obtained from a differential system

\( \begin{align} \frac{dy(t)}{dt}+\frac{y(t)}{4}&=x(t) \end{align} \)

You first need to rearrange this to be an equation for \(\frac{dy(t)}{dt}\) as a function of \(t\), \(y\), and constants. This equation also has an independent \(x(t)\) so that will be a part of the function as well:

\( \begin{align} \frac{dy(t)}{dt}&=x(t)-\frac{y(t)}{4}\\ \end{align} \)

Now if you further assume that:

\( \begin{align} x(t)&=\cos(3t)\\ y(0)&=-1 \end{align}\)

you can numerically solve this differential equation (and plot the results for both \(x(t)\) and \(y(t)\)) with:

# %% Imports
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# %% Define independent function and derivative function
def x(t):
    return np.cos(3 * t)

def f(t, y, c):
    dydt = x(t) - y / 4
    return dydt

# %% Define time spans, initial values, and constants
tspan = np.linspace(0, 15, 1000)
yinit = [-1]
c = []

# %% Solve differential equation
sol = solve_ivp(lambda t, y: f(t, y, c),
                [tspan[0], tspan[-1]], yinit, t_eval=tspan, rtol = 1e-5)

# %% Plot independent and dependent variable
# Note that sol.y[0] is needed to extract a 1-D array
plt.figure(1)
plt.clf()
fig, ax = plt.subplots(num=1)
ax.plot(sol.t, x(sol.t), 'k-', label='Input')
ax.plot(sol.t, sol.y[0], 'k--', label='Output')
ax.legend(loc='best')

As an aside, to get the same graph in Maple, you could use:

restart;

deqn := diff(y(t), t) = x(t)-(1/4)*y(t);

Src := x(t) = cos(3*t);

MySoln := dsolve({subs(Src, deqn), y(0) = 3}, [y(t)]);

plot(subs(Src, MySoln, [x(t), y(t)]), t = 0 .. 15);

More Examples

There are several more examples at Python:Ordinary_Differential_Equations/Examples -- these are all based on starter code at Python:Ordinary_Differential_Equations/Templates

Notes / Troubleshooting

• The initial \(y\) value is assume to be at the start of the time span t_span, which means not neessarily time 0.
• The initial values need to be in "an array-like" thing; either a list or an array. This is true even if there is a single initial condition. If you are getting an error:

`y0` must be 1-dimensional.

that means you need to enclose the initial guesses in brackets.

Questions

Post your questions by editing the discussion page of this article. Edit the page, then scroll to the bottom and add a question by putting in the characters *{{Q}}, followed by your question and finally your signature (with four tildes, i.e. ~~~~). Using the {{Q}} will automatically put the page in the category of pages with questions - other editors hoping to help out can then go to that category page to see where the questions are. See the page for Template:Q for details and examples.

External Links

SciPy page for solve_ivp