Python:Ordinary Differential Equations

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 the first function in concert with SciPy's ODE solvers to calculate solutions over a specified time range assuming given initial conditions.

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}&=x(t)-\frac{y(t)}{4}\\ x(t)&=\cos(3t)\\ y(0)&=5 \end{align}\)

you can write that code as:

#%% 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):
    dydt = x(t)-y/4
    return dydt
    
#%% Define time spans and initial value
tspan = np.linspace(0, 15, 1000)
yinit = [5]

#%% Solve differential equation
sol = solve_ivp(f, [tspan[0], tspan[-1]], yinit, t_eval=tspan)

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

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) = 5}, [y(t)]);

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

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