IVP systems - Fundamentals of Numerical Computation

Few applications involve an initial-value problem with just a single dependent variable. Usually there are multiple unknowns and a system of equations to define them.

Example 6.3.1

Variations of the following model are commonly seen in ecology:

\begin{split} \frac{d y}{d t} &= y(1-\alpha y) - \frac{yz}{1+\beta y}, \\ \frac{d z}{d t} &= -z + \frac{yz}{1+\beta y}, \end{split}

(6.3.1)

where α and β are positive constants. This model is a system of two differential equations for the unknown functions $y(t)$ , which could represent a prey species or susceptible host, and $z(t)$ , which could represent a predator species or infected population. We refer to this as a predator–prey model. Both of the equations involve both of the unknowns, with no clear way to separate them.

We can pack the two dependent variables $y$ and $z$ into a vector-valued function of time, $\mathbf{u}(t)$ , writing

\begin{split} u_1'(t) &= f_1(t,\mathbf{u}) = u_1(1-au_1) - \frac{u_1 u_2}{1+bu_1},\\ u_2'(t) &= f_2(t,\mathbf{u}) = -u_2 + \frac{u_1 u_2}{1+bu_1}, \end{split}

(6.3.2)

and identifying $u_1=y$ , $u_2=z$ .

We now upgrade our IVP definition, Definition 6.1.1.

We use the terms IVP system and vector-valued IVP interchangeably; a system of scalar IVPs can be put into the form of (6.3.3) by appropriate definitions of $\mathbf{u}$ and $\mathbf{f}$ , as shown in Example 6.3.1.

6.3.1Numerical solutions¶

The generalization of any scalar IVP solver to handle systems is straightforward. Consider Euler’s method, which in system form becomes

\begin{split} \mathbf{u}_{i+1} &= \mathbf{u}_i + h\,\mathbf{f}(t_i,\mathbf{u}_i), \qquad i=0,\ldots,n-1. \end{split}

(6.3.4)

The vector difference equation (6.3.4) is just Euler’s formula applied simultaneously to each component of the ODE system. Because operations such as addition and multiplication translate easily from scalars to vectors, Function 6.2.2 that we wrote for scalar IVPs works for systems as well. Practically speaking, the only changes that must be made are that the initial condition and the ODE function have to be coded to use vectors.

Example 6.3.2 (Predator-prey model)

Julia

MATLAB

Python

Example 6.3.2

We encode the predator–prey equations via a function.

function predprey(u, p, t)
    α, β = p      # rename parameters for convenience
    y, z = u      # rename solution components
    s = (y * z) / (1 + β * y)     # appears in both equations
    return [y * (1 - α * y) - s, -z + s]
end;

As before, the ODE function must accept three inputs, u, p, and t, even though in this case there is no explicit dependence on t. The second input is used to pass parameters that don’t change throughout a single instance of the problem.

To specify the IVP we must also provide the initial condition, which is a 2-vector here, and the interval for the independent variable.

using OrdinaryDiffEq
u₀ = [1, 0.01]
tspan = (0.0, 60.0)
α, β = 0.1, 0.25
ivp = ODEProblem(predprey, u₀, tspan, [α, β])

ODEProblem with uType Vector{Float64} and tType Float64. In-place: false
timespan: (0.0, 60.0)
u0: 2-element Vector{Float64}:
 1.0
 0.01

You can use any DifferentialEquations solver on the IVP system.

using Plots
sol = solve(ivp, Tsit5());
plot(sol, label=["prey" "predator"],
    title="Predator-prey solution")

We can find the discrete values used to compute the interpolated solution. The sol.u value is a vector of vectors.

t, u = sol.t, sol.u    # extract times and solution values
@show size(u);
@show t[20];
@show u[20];

size(u) = (100,)


t[20] = 8.386698083051241
u[20] = [0.027739191695661334, 0.699418322938034]

We can also use Function 6.2.2 to find the solution.

t, u = FNC.euler(ivp, 1200);

The solution u is a vector of [prey,predator] 2-vectors for each of the discrete times in t. Manipulating the vector-of-vectors output can be a little tricky. Here, we convert it to an $n\times 2$ matrix. Each column is one component, while each row is a single value of $t$ .

u = [u[j] for u in u, j in 1:2]
plot!(t[1:3:end], u[1:3:end, :];
    l=(1, :black),  m=2,
    label=["Euler prey" "Euler predator"])

Notice above that the accuracy of the Euler solution deteriorates rapidly.

When there are just two components, it’s common to plot the solution in the phase plane, i.e., with $u_1$ and $u_2$ along the axes and time as a parameterization of the curve.

plot(sol, idxs=(1, 2),
    title="Predator-prey in the phase plane",
    xlabel=L"y",  ylabel=L"z")

From this plot we can deduce that the solution approaches a periodic one, which in the phase plane is represented by a closed loop.

Example 6.3.2

We encode the predator–prey equations via a function, defined here externally.

f63_predprey.m

function du_dt = predprey(t, u, p)
    alpha = p(1);  beta = p(2);
    y = u(1);      z = u(2);
    s = (y * z) / (1 + beta * y);  % appears in both equations
    du_dt = [ y * (1 - alpha * y) - s;  -z + s ];
end

The values of alpha and beta are parameters that influence the solution of the IVP. We use the Parameters field of the IVP object to define them for the solver, which in turn passes them as the third argument into our ODE function.

u0 = [1; 0.01];    % column vector
p = [0.1, 0.25];
ivp = ode;
ivp.ODEFcn = @f63_predprey;
ivp.InitialValue = u0;
ivp.Parameters = p;
sol = solve(ivp, 0, 60);
size(sol.Solution)

Each column of the Solution field is the solution vector $\mathbf{u}$ at a particular time; each row is a component of $\mathbf{u}$ over all time.

clf
plot(sol.Time, sol.Solution)
xlabel("t")
ylabel("u(t)")
title('Predator-prey solution')
legend('prey', 'predator');

We can also use Function 6.2.2 to find the solution.

[t, u] = eulerivp(ivp, 0, 60, 1200);

hold on
plot(t, u, '.')

Notice above that the accuracy of the Euler solution deteriorates rapidly.

When there are just two components, it’s common to plot the solution in the phase plane, i.e., with $u_1$ and $u_2$ along the axes and time as a parameterization of the curve.

clf
plot(u(1, :), u(2, :))
title("Predator-prey in the phase plane")
xlabel("y")
ylabel(("z"));

From this plot we can deduce that the solution approaches a periodic one, which in the phase plane is represented by a closed loop.

Example 6.3.2

We encode the predator–prey equations via a function.

def predprey(t, u):
    y, z = u                        # rename for convenience
    s = (y * z) / (1 + beta * y)    # appears in both equations
    return array([y * (1 - alpha * y) - s, -z + s])

To specify the IVP we must also provide the initial condition, which is a 2-vector here, and the interval for the independent variable. These are given in the call to solve_ivp.

from scipy.integrate import solve_ivp
u0 = array([1, 0.01])
tspan = [0.0, 80.0]
alpha, beta = 0.1, 0.25
sol = solve_ivp(predprey, tspan, u0, dense_output=True)
print(f"solved with {sol.y.shape[1]} time steps")

solved with 102 time steps

As in scalar problems, the solution object has fields t and y that contain the values of the independent and dependent variables, respectively. Each row of y represents one component of the solution at every time step, and each column of y is the entire solution vector at one time step. Since we used dense_output=True, there is also a method sol that can be used to evaluate the solution at any time.

t = linspace(0, 80, 1200)
u = vstack([sol.sol(t[i]) for i in range(t.size)]).T    # same shape as sol.y
fig, ax = subplots()
ax.plot(t, u[0, :], label="prey")
ax.plot(t, u[1, :], label="predator")
xlabel("$t$"), ylabel("population")
title(("Predator-prey solution"));

We can also use Function 6.2.2 to find the solution.

t_E, u_E = FNC.euler(predprey, tspan, u0, 800)
ax.scatter(t_E, u_E[0, :], label="prey (Euler)", s=1)
ax.scatter(t_E, u_E[1, :], label="predator (Euler)", s=2)
ax.legend()
fig

You can see above that the Euler solution is not very accurate. When the solution has two components, it’s common to plot the it in the phase plane, i.e., with $u_1$ and $u_2$ along the axes and time as a parameterization of the curve.

plot(u[0, :], u[1, :])
xlabel("prey"), ylabel("predator")
title(("Predator-prey phase plane"));

From this plot we can see that the solution approaches a periodic one, which in the phase plane is represented by a closed path.

In the rest of this chapter we present methods as though they are for scalar equations, but their application to systems is taken for granted. The generalization of error analysis can be more complicated, but our statements about order of accuracy and other properties are true for systems as well as scalars. The codes are all written to accept systems.

6.3.2Transformation of high-order systems¶

Fortunately, the ability to solve first-order ODE systems implies the ability to solve systems of higher differential order, too. The reason is that there is a systematic way to turn a higher-order problem into a first-order one of higher dimension.

Example 6.3.4

Two identical pendulums suspended from the same rod and swinging in parallel planes can be modeled as the second-order system

\begin{split} \theta_1''(t) +\gamma \theta_1' + \frac{g}{L} \sin \theta_1 + k(\theta_1-\theta_2) &= 0,\\ \theta_2''(t) +\gamma \theta_2' + \frac{g}{L} \sin \theta_2 + k(\theta_2-\theta_1) &= 0, \end{split}

(6.3.8)

where $\theta_1$ and $\theta_2$ are angles made by the two pendulums, $L$ is the length of each pendulum, γ is a frictional parameter, and $k$ is a parameter describing a torque produced by the rod when it is twisted. We can convert this problem into a first-order system using the substitutions

u_1 = \theta_1, \quad u_2 = \theta_2, \quad u_3 = \theta_1', \quad u_4 = \theta_2'.

(6.3.9)

With these definitions the system becomes

\begin{split} u_1' &= u_3, \\ u_2' &= u_4, \\ u_3' &= -\gamma u_3 - \frac{g}{L}\sin u_1 + k(u_2-u_1), \\ u_4' &= -\gamma u_4 - \frac{g}{L}\sin u_2 + k(u_1-u_2), \end{split}

(6.3.10)

which is a first-order system in four dimensions. To complete the description of the problem, you need to specify values for $\theta_1(0)$ , $\theta_1'(0)$ , $\theta_2(0)$ , and $\theta_2'(0)$ .

The trick illustrated in the preceding examples is always available. Suppose $y$ is a scalar dependent variable in the system. You should introduce a component of $\mathbf{u}$ for $y$ , $y'$ , etc., up to but not including the highest derivative appearing anywhere for $y$ . This is done for each scalar variable in the original system. There should be one component of $\mathbf{u}$ for each scalar initial condition given. Many equations for the first-order system then come from the trivial relationships among all the lower derivatives. The remaining equations for the system come from the original, high-order equations. In the end, there must be as many scalar component equations as unknown first-order variables.

Example 6.3.5 (Coupled pendulums)

Julia

MATLAB

Python

Example 6.3.5

Let’s implement the coupled pendulums from Example 6.3.4. The pendulums will be pulled in opposite directions and then released together from rest.

function couple(u, p, t)
    γ, L, k = p
    g = 9.8
    udot = similar(u)
    udot[1:2] .= u[3:4]
    udot[3] = -γ * u[3] - (g / L) * sin(u[1]) + k * (u[2] - u[1])
    udot[4] = -γ * u[4] - (g / L) * sin(u[2]) + k * (u[1] - u[2])
    return udot
end

u₀ = [1.25, -0.5, 0, 0]
tspan = (0.0, 50.0);

First we check the behavior of the system when the pendulums are uncoupled, i.e., when $k=0$ .

γ, L, k = 0, 0.5, 0
ivp = ODEProblem(couple, u₀, tspan, [γ, L, k])
sol = solve(ivp, Tsit5())
plot(sol, idxs=[1, 2], 
    label=[L"\theta_1" L"\theta_2"],
    xlims=[20, 50], 
    title="Uncoupled pendulums")

You can see that the pendulums swing independently. Because the model is nonlinear and the initial angles are not small, they have slightly different periods of oscillation, and they go in and out of phase.

With coupling activated, a different behavior is seen.

k = 1
ivp = ODEProblem(couple, u₀, tspan, [γ, L, k])
sol = solve(ivp, Tsit5())
plot(sol, idxs=[1, 2], 
    label=[L"\theta_1" L"\theta_2"],
    xlims=[20, 50], 
    title="Coupled pendulums")

The coupling makes the pendulums swap energy back and forth.

Example 6.3.5

Let’s implement the coupled pendulums from Example 6.3.4. The pendulums will be pulled in opposite directions and then released together from rest.

f63_pendulums.m

function udot = pendulums(t, u, p)
    gamma = p(1);  L = p(2);  k = p(3);
    g = 9.8;
    udot = zeros(4, 1);
    udot(1:2) = u(3:4);
    udot(3) = -gamma * u(3) - (g / L) * sin(u(1)) + k * (u(2) - u(1));
    udot(4) = -gamma * u(4) - (g / L) * sin(u(2)) + k * (u(1) - u(2));
end

u0 = [1.25; -0.5; 0; 0];
a = 0; b = 50;

First we check the behavior of the system when the pendulums are uncoupled, i.e., when $k=0$ .

params =[0.01, 0.5, 0];    % gamma, L, k
ivp = ode(ODEFcn=@f63_pendulums, InitialValue=u0, Parameters=params);
theta = solutionFcn(ivp, a, b, OutputVariables = 1:2);
t = linspace(a, b, 1001);
clf, plot(t, theta(t))
xlabel("t");  ylabel("angle")
title("Uncoupled pendulums")
legend("\theta_1", "\theta_2");

With coupling activated, a different behavior is seen.

params(3) = 1;
ivp = ode(ODEFcn=@f63_pendulums, InitialValue=u0, Parameters=params);
theta = solutionFcn(ivp, a, b, OutputVariables = 1:2);
clf, plot(t, theta(t))
xlabel("t");  ylabel("angle")
title("Coupled pendulums")
legend("\theta_1", "\theta_2");

The coupling makes the pendulums swap energy back and forth.

Example 6.3.5

Let’s implement the coupled pendulums from Example 6.3.4. The pendulums will be pulled in opposite directions and then released together from rest.

def couple(t, u, params):
    gamma, L, k = params
    g = 9.8
    udot = copy(u)
    udot[:2] = u[2:4]
    udot[2] = -gamma * u[2] - (g / L) * sin(u[0]) + k * (u[1] - u[0])
    udot[3] = -gamma * u[3] - (g / L) * sin(u[1]) + k * (u[0] - u[1])
    return udot

u0 = array([1.25, -0.5, 0, 0])
tspan = [0.0, 50.0]

First we check the behavior of the system when the pendulums are uncoupled, i.e., when $k=0$ .

gamma, L, k = 0.01, 0.5, 0.0
du_dt = lambda t, u: couple(t, u, (gamma, L, k))
sol = solve_ivp(du_dt, tspan, u0, t_eval=linspace(0, 50, 1000))
plot(sol.t, sol.y[:2, :].T)    # first two components of solution
xlabel("t"), ylabel("angle")
title("Uncoupled pendulums");

With coupling activated, a different behavior is seen.

k = 0.75    # changes the value in the du_dt closure
sol = solve_ivp(du_dt, tspan, u0, t_eval=linspace(0, 50, 1000))
plot(sol.t, sol.y[:2, :].T)
xlabel("t"), ylabel("angle")
title("Coupled pendulums");

The coupling makes the pendulums swap energy back and forth.

6.3.3Exercises¶

✍ Rewrite the given higher order problems as first-order systems.
(a) $y'''-3y''+3 y' -y = t, \: y(0) = 1, \: y'(0) = 2, \: y''(0) = 3$
(b) $y'' + 4 (x^2-1)y' + y = 0, \: y(0) = 2, \: y'(0) = -1$
(c) For a given constant $a$ ,
$\begin{split} x'' + \frac{a x}{(x^2+y^2)^{3/2}} &= 0,\\ y'' + \frac{a y}{(x^2+y^2)^{3/2}} &= 0, \end{split}$
(6.3.11)
with initial values $x(0) = 1$ , $x'(0)=y(0) = 0$ , $y'(0)=3$
(d) $y^{(4)} -y = e^{-t}, \: y(0) = 0, \: y'(0) = 0, \: y''(0) = 1,\: y'''(0) = 0$
(e) $y'''-y''+y'-y = t, \: y(0) = 1, \: y'(0) = 2, \: y''(0) = 3$
✍ Write the given IVP as a system. Then do two steps of Euler’s method by hand (perhaps with a calculator) with the indicated step size $h$ . Using the given exact solution, compute the error after the second step.
(a) $y''+ 4y = 4t, \: y(0) = 1,\: y'(0) = 1; \: \hat{y}(t) = t+\cos (2t),\: h=0.1$
(b) $y''- 4y = 4t, \: y(0) = 2,\: y'(0) = -1; \: \hat{y}(t) = e^{2t} + e^{-2t}-t,\: h=0.1$
(c) $2 x^2 y'' +3xy' - y = 0, \: y(2) = 1, \: y'(2) = -1/2, \: \hat{y}(x) = 2/x, h = 1/8$
(d) $2 x^2 y'' +3xy' - y = 0,\: y(1) = 4, \: y'(1) = -1, \: \hat{y}(x) = 2(x^{1/2} + x^{-1}), h=1/4$
⌨ Solve the following IVPs using Function 6.2.2 using $n=1000$ steps. Plot the solution and its first derivative together on one plot, and plot the error in each component as functions of time on another.
(a) $y''+ 4y = 4t, \: 0< t< 2\pi, \: y(0) = 1,\: y'(0) = 1; \: \hat{y}(t) = t+\cos (2t)$
(b) $y''+ 9y = \sin(2t), \: 0< t< 2\pi, \: y(0) = 2,\: y'(0) = 1$ ; $\quad \hat{y}(t) = (1/5) \sin(3t) + 2 \cos (3t)+ (1/5) \sin (2t)$
(c) $y''- 4y = 4t \: 0< t< 1.5, \: y(0) = 2,\: y'(0) = -1; \: \hat{y}(t) = e^{2t} + e^{-2t}-t$
(d) $y''+ 4y'+ 4y = t, \: 0< t< 4, \: y(0) = 1,\: y'(0) = 3/4; \: \hat{y}(t) = (3t+5/4)e^{-2t} + (t-1)/4$
(e) $x^2 y'' +5xy' + 4y = 0,\: 1<x<e^2, \: y(1) = 0, \: y'(1) = 2, \: \hat{y}(x) = (2/x^2) \ln x$
(f) $x^2 y'' +5xy' + 4y = 0,\: 1<x<e^2, \: y(1) = 1, \: y'(1) = -1, \: \hat{y}(x) = x^{-2}( 1 + \ln x)$
(g) $2 x^2 y'' +3xy' - y = 0,\: 2<x<20, \: y(2) = 1, \: y'(2) = -1/2, \: \hat{y}(x) = 2/x$
(h) $2 x^2 y'' +3xy' - y = 0,\: 1<x<16, \: y(1) = 4, \: y'(1) = -1, \: \hat{y}(x) = 2(x^{1/2} + x^{-1})$
(i) $x^2 y'' -xy' + 2y = 0,\: 1<x<e^{\pi}, \: y(1) = 3, \: y'(1) = 4$ ; $\quad \hat{y}(x) = x \left[ 3 \cos \left( \ln x \right)+\sin \left( \ln x \right) \right]$
(j) $x^2 y'' + 3xy' + 4y = 0,\: e^{\pi/12} < x < e^{\pi}, \: y(e^{\pi/12}) = 0, \: y'(e^{\pi/12}) = -6$ ; $\quad \hat{y}(x) = x^{-1} \left[ 3 \cos \left( 3 \ln x \right)+\sin \left( 3 \ln x \right) \right]$

⌨ A disease that is endemic to a population can be modeled by tracking the fraction of the population that is susceptible to infection, $v(t)$ , and the fraction that is infectious, $w(t)$ . (The rest of the population is considered to be recovered and immune.) A typical model is the SIR model (see Britton (2003))
$\frac{dv}{dt} = 0.2(1-v) - 3vw, \qquad \frac{dw}{dt} = (3v-1)w.$
(6.3.12)
Starting with $v(0) = 0.95$ and $w(0) = 0.05$ , use solve to find the long-term steady values of $v(t)$ and $w(t)$ . Plot both components of the solution as functions of time.
⌨ In each case below, use solve to solve the given ODE for $0\le t \le 10$ with the given initial conditions. Plot the results together as curves in the phase plane (that is, with $x$ and $y$ as the axes of the plot), using aspect_ratio=1 in the plot command.
(a)
$\begin{split} x'(t) & = - 4y + x(1-x^2-y^2),\\ y'(t) & = 4x + y(1-x^2-y^2), \end{split}$
(6.3.13)
with $[x(0),y(0)]=[0.1,0]$ and $[x(0),y(0)]=[0,1.9]$ .
(b)
$\begin{split} x'(t) & = - 4y - \tfrac{1}{4}x(1-x^2-y^2)(4-x^2-y^2),\\ y'(t) & = 4x - \tfrac{1}{4}y(1-x^2-y^2)(4-x^2-y^2), \end{split}$
(6.3.14)
with $[x(0),y(0)]=[0.95,0]$ , $[0,1.05]$ , and $[-2.5,0]$ .

⌨ The FitzHugh–Nagumo equations are a simple model of the repeated firing of a neuron. They are given by
$\begin{split} \frac{d v_1}{dt} &= - v_1(v_1-1)(v_1-a) - v_2 + I, \\ \frac{d v_2}{dt} &= \epsilon ( v_1 - \gamma v_2). \end{split}$
(6.3.15)
Assume $v_1(0) = 0.5$ , $v_2(0) = 0.1$ , $a = 0.1$ , $\epsilon = 0.008$ , $\gamma = 1$ . For each value of $I$ below, find and plot the solution using solve for $0\le t \le 600$ . The solutions are highly sensitive to $I$ , and you need to change the requested absolute and relative error tolerances to 10^-9. In each case the solution quickly approaches a periodic oscillation.
(a) $I = 0.05527,\quad$ (b) $I = 0.05683,\quad$ (c) $I = 0.0568385,\quad$ (d) $I = 0.05740$ .
This exploration was carried out by Baer and Erneux Baer & Erneux (1986).

References¶

Britton, N. F. (2003). Essential Mathematical Biology. Springer.
Baer, S. M., & Erneux, T. (1986). Singular Hopf Bifurcation to Relaxation Oscillations. SIAM Journal on Applied Mathematics, 46(5), 721–739. 10.1137/0146047