How to do coupled and 2nd order ODEs

13.1 ODEs: Coupled

Consider these two differential equations:

$\begin{array}{rcl} \displaystyle\frac{dy_0}{dt}&=&y_1\\ &&\\ \displaystyle\frac{dy_1}{dt}&=&-y_0-\displaystyle\frac{1}{10}y_1 \end{array}$

Since the evolution in time of $y_0(t)$ depends on $y_1(t)$ and vice versa, the equations constitute a set of coupled differential equations.

Given some initial value condition, e.g.

$\begin{array}{rcl} y_0(t=0)&=&1\\ &&\\ y_1(t=0)&=&-1 \end{array}$

the problem of finding $y_0(t)$ and $y_1(t)$ until some later time, e.g. $t_\mathrm{final}=20$ , may be solved with solve_ivp like this:

from scipy.integrate import solve_ivp

tinit = 0
tfinal = 20
trange = [tinit,tfinal]
yinit = [1,-1]

dydt = lambda t, y: [y[1],-y[0]-1/10*y[1]]
mysol = solve_ivp(dydt, trange, yinit)
ts = mysol.t
y0s = mysol.y[0]
y1s = mysol.y[1]

Compared to the situation with a single ODE, a number of things has happened:

the function dydt now assumes that its 2nd argument is a list with two elements, $y_0$ and $y_1$ .
the function dydt returns a list of the right-hand sides of both differential equations, i.e. it returns $\Bigl[\displaystyle\frac{dy_0}{dt},\frac{dy_1}{dt}\Bigr]$ .
the initial value for the dependent variable, yinit, is now a list $\Bigl[y_0,y_1\Bigr]$
two lists of solutions (for $y_0$ and $y_1$ , respectively) are extracted from the result of solve_ivp with the mysol.y[0] and mysol.y[1] statements.

When plotted with these commands:

import matplotlib.pyplot as plt
plt.rc('font', size=16)

fig, ax = plt.subplots()
ax.plot(ts,y0s,linestyle='dashed',marker='o')
ax.plot(ts,y1s,linestyle='dashed',marker='x')
ax.set_xlabel('$t$')
ax.set_ylabel('$y$')
ax.legend(['$y_0$','$y_1$'])

the solution looks like this

Figure 1

13.2 Different naming of variables

Consider if we had encounted the exact same set of differential equations, only with another naming of the dependent variables, e.g.:

$\begin{array}{rcl} \displaystyle\frac{dr}{dt}&=&w\\ &&\\ \displaystyle\frac{dw}{dt}&=&-r-\displaystyle\frac{1}{10}w \end{array}$

with the initial value condition: $r=1$ and $w=-1$ for $t=0$ . Well, obviously, the solution would be the same, and the above Python program would do the job. However, it would be neat the write the program with a naming convention for the variables that is closer to the problem as it is formulated. That could for instance be achieved by writing:

tinit = 0
tfinal = 20
trange = [tinit,tfinal]
rinit = 1
winit = -1
yinit = [rinit,winit]

def dydt(t,y):
    r = y[0]
    w = y[1]
    drdt = w
    dwdt = -r - 1/10 * w
    return [drdt,dwdt]
    
mysol = solve_ivp(dydt, trange, yinit)
ts = mysol.t
rs = mysol.y[0]
ws = mysol.y[1]

which provides the exact same solution:

Figure 2

Note how the function dydt evaluating the right-hand side of the differential equations has been defined at length, rather than in the compact form of a lambda-function. This in particular allows for the section,

    r = y[0]
    w = y[1]

where the original naming of variables is recoved, so that when the right-hand sides of the differential equations are evaluated, they become readable:

    drdt = w
    dwdt = -r - 1/10 * w

and hence immediately recognizable compared to the problem formulation.

13.3 ODE: 2nd order

A 2nd order differential equation with one dependent variable generally looks like this:

$\frac{d^2x}{dt^2}=f(x,\frac{dx}{dt}),$

where $f(x,\displaystyle\frac{dx}{dt})$ may be anything. Consider as an example that we want to solve:

$\frac{d^2x}{dt^2}=-x-\frac{1}{10}\frac{dx}{dt}$

and that the initial value condition is specified as $x(t=0)=1$ and $\displaystyle\left.\frac{dx}{dt}\right|_{t=0}=-1$ . Introducing a "helper"-variable,

$v=\displaystyle\frac{dx}{dt},$

we may rewrite the 2nd order differential equation as two 1st order differential equations without loss of generality:

$\begin{array}{rcl} \displaystyle\frac{dx}{dt}&=&v\\ \\ \displaystyle\frac{dv}{dt}&=&-x-\displaystyle\frac{1}{10}v \end{array}$

whereby the initial value condition becomes $x(t=0)=1$ and $v(t=0)=-1$ .

As a consequence of the rewriting, the 2nd order differential equation has turned into a set of coupled 1st order differential equations, that looks exactly like some we have already solved with solve_ivp. Adopting the same code, but using the present variable names, it may look like this:

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

tinit = 0
tfinal = 12
trange = [tinit,tfinal]
x0 = 1
v0 = -1
yinit = [x0,v0]

def dydt(t,y):
    x = y[0]
    v = y[1]
    dxdt = v
    dvdt = -x - 1/10 * v
    return [dxdt,dvdt]

mysol = solve_ivp(dydt, trange, yinit)
ts = mysol.t
xs = mysol.y[0]
vs = mysol.y[1]

plt.rc('font', size=16)
fig,ax = plt.subplots(1,1)
ax.plot(ts,xs,'o',linestyle='--')
ax.plot(ts,vs,'x',linestyle='--')
ax.grid()
ax.set_xlim([-1,13])
ax.set_xlabel('$t$')
ax.set_ylabel('$x,v$')
ax.legend(['$x$',r'$v=\frac{dx}{dt}$'])

whereby we have solved a 2nd order differential equation. The plot looks like this:

Figure 3