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.1 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
Non-trivial mass matrix: 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.1 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 solution in the phase plane",
    xlabel=L"y(t)",  ylabel=L"z(t)")

As time progresses, the point in the phase plane spirals inward toward a limiting closed loop called a limit cycle representing a periodic solution:

Source

prey, predator = [], []
anim = @animate for t in range(0, 60, 801)
    y, z = sol(t)
    plot(prey, predator;
        title="Predator-prey solution in the phase plane",
        xlabel=L"y(t)",  ylabel=L"z(t)", legend=false)
    push!(prey, y)
    push!(predator, z)
    scatter!([y], [z], m=(5, :red))
    xlims!(0, 9)
    ylims!(0, 6)
end
mp4(anim, "figures/predator-prey.mp4")

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, so the function signature above includes a third argument for them. But the solvers we use expect to receive a function of $t$ and $\mathbf{u}$ only. To deal with this, we create a function handle f that captures the parameter values in its definition.

params = [0.1, 0.25];
f = @(t, u) f63_predprey(t, u, params);

We use a new syntax here for calling ode45. By giving it a vector as the tspan argument, we are asking it to return the solution at the times in that vector.

u0 = [1; 0.01];    % column vector
a = 0; b = 60;     % time interval

[t, u] = ode45(f, linspace(a, b, 1001), u0);
size(u)

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

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

We can also use Function 6.2.1 to find the solution (although it does not allow us to choose different time nodes for the output).

[t_e, u_e] = eulersys(f, [a, b], u0, 1200);

hold on
plot(t_e, u_e, '.', displayname='Euler')

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 solution in the phase plane")
xlabel("prey, y(t)")
ylabel("predator, z(t)");

As time progresses, the point in the phase plane spirals inward toward a limiting closed loop called a limit cycle representing a periodic solution:

Source

axis(axis)
cla,  hold on
grid on
curve_ = plot(u(1, 1), u(1, 2))
head_ = plot(u(1, 1), u(1, 2), 'o', 'markersize', 10, 'markerfacecolor', 'r');
title('Predator–prey solution in the phase plane') 
xlabel("prey, y(t)")
ylabel("predator, z(t)");
vid = VideoWriter("figures/predator-prey.mp4","MPEG-4");
vid.Quality = 85;
open(vid)
for frame = 1:2:length(t)
    curve_.XData = u(1:frame, 1);
    curve_.YData = u(1:frame, 2);
    head_.XData = u(frame, 1);
    head_.YData = u(frame, 2);
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)

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, 60.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 76 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, 60, 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.1 to find the solution.

t_E, u_E = FNC.euler(predprey, tspan, u0, 1001)
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

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

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, $y(t)$"), ylabel("predator, $z(xt)$")
title("Predator-prey solution in the phase plane");

As time progresses, the point in the phase plane spirals inward toward a limiting closed loop called a limit cycle representing a periodic solution:

Source

from matplotlib import animation
fig, ax = subplots()
curve_ = ax.plot(u[0, 0], u[1, 0])[0]
head_ = ax.plot(u[0, 0], u[1, 0], 'ro', markersize=6)[0]

ax.set_xlabel("prey, $y(t)$")
ax.set_ylabel("predator, $z(xt)$")
ax.set_xlim(0, 9)
ax.set_ylim(0, 6)
ax.set_title("Predator-prey solution in the phase plane")

predator = []
prey = []
def snapshot(t):
    y, z = sol.sol(t)
    prey.append(y)
    predator.append(z)
    head_.set_data([y], [z])
    curve_.set_data(prey, predator)

anim = animation.FuncAnimation(
    fig, snapshot, frames=linspace(0, 60, 801)
    )
anim.save("figures/predator-prey.mp4", fps=24)
close()

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.3, 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:

Source

anim = @animate for t in range(0, 24, 251)
    θ₁, θ₂, _ = sol(t)
    plot([0, sin(θ₁)], [0, -cos(θ₁)], l=4;
        layout=(1, 2), aspect_ratio=1, grid=false,
        xaxis=((-1.1, 1.1), false), yaxis=((-1.1, 0.1), false), legend=false)
    plot!([0, sin(θ₂)], [0, -cos(θ₂)], l=4;
        subplot=2, aspect_ratio=1, grid=false,
        xaxis=((-1.1, 1.1), false), yaxis=((-1.1, 0.1), false), legend=false)
    scatter!([sin(θ₁)], [-cos(θ₁)], m=(5, :black), subplot=1)
    scatter!([sin(θ₂)], [-cos(θ₂)], m=(5, :black), subplot=2)
    annotate!(-0.95, 0.25, Plots.text(@sprintf("t = %.1f", t), :left, 11); subplot=1)
end
mp4(anim, "figures/pendulums-weak.mp4")

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:

Source

anim = @animate for t in range(0, 50, 601)
    θ₁, θ₂, _ = sol(t)
    plot([0, sin(θ₁)], [0, -cos(θ₁)], l=4;
        layout=(1, 2), aspect_ratio=1, grid=false,
        xaxis=((-1.1, 1.1), false), yaxis=((-1.1, 0.1), false), legend=false)
    plot!([0, sin(θ₂)], [0, -cos(θ₂)], l=4;
        subplot=2, aspect_ratio=1, grid=false,
        xaxis=((-1.1, 1.1), false), yaxis=((-1.1, 0.1), false), legend=false)
    scatter!([sin(θ₁)], [-cos(θ₁)], m=(5, :black), subplot=1)
    scatter!([sin(θ₂)], [-cos(θ₂)], m=(5, :black), subplot=2)
    annotate!(-0.95, 0.1, Plots.text(@sprintf("t = %.1f", t), :left, 11); subplot=1)
end
mp4(anim, "figures/pendulums-strong.mp4")

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.3; 0; 0];
tspan = [0, 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
f = @(t, u) f63_pendulums(t, u, params);
sol = ode45(f, tspan, u0);
theta = @(t) deval(sol, t, 1:2)';
t = linspace(tspan(1), tspan(2), 901);
clf, plot(t, theta(t))
xlabel("t");  ylabel("angle")
title("Uncoupled pendulums")
legend("\theta_1", "\theta_2");

You can see that the pendulums swing independently:

Source

clf
rod_ = [];  bob_ = [];
for k = 1:2
    subplot(1, 2, k)
    axis equal, axis([-1.1, 1.1, -1.1, 1.1])
    hold on, axis off
    rod_ = [rod_; plot(0, 0, 'linewidth', 3)];
    bob_= [bob_; plot(0, 0, 'ko', 'markersize', 20, 'markerfacecolor', 'k')];
end
text_ = text(-0.9, 0.85, "t = 0", 'fontsize', 20);
vid = VideoWriter("figures/pendulums-weak.mp4","MPEG-4");
vid.FrameRate = 24;
open(vid)
for frame = 1:401
    q = theta(t(frame));
    for k = 1:2
        rod_(k).XData = [0, sin(q(k))];
        rod_(k).YData = [0, -cos(q(k))];
        bob_(k).XData = sin(q(k));
        bob_(k).YData = -cos(q(k));
    end
    text_.String = sprintf("t = %.2f", t(frame));
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)

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.

params(3) = 1;
f = @(t, u) f63_pendulums(t, u, params);
sol = ode45(f, tspan, u0);
theta = @(t) deval(sol, t, 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:

Source

clf
rod_ = [];  bob_ = [];
for k = 1:2
    subplot(1, 2, k)
    axis equal, axis([-1.1, 1.1, -1.1, 1.1])
    hold on, axis off
    rod_ = [rod_; plot(0, 0, 'linewidth', 3)];
    bob_= [bob_; plot(0, 0, 'ko', 'markersize', 20, 'markerfacecolor', 'k')];
end
text_ = text(-0.9, 0.85, "t = 0", 'fontsize', 20);
vid = VideoWriter("figures/pendulums-strong.mp4","MPEG-4");
vid.FrameRate = 24;
open(vid)
for frame = 1:length(t)
    q = theta(t(frame));
    for k = 1:2
        rod_(k).XData = [0, sin(q(k))];
        rod_(k).YData = [0, -cos(q(k))];
        bob_(k).XData = sin(q(k));
        bob_(k).YData = -cos(q(k));
    end
    text_.String = sprintf("t = %.2f", t(frame));
    writeVideo(vid, frame2im(getframe(gcf)));
end
close(vid)

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, 1001))
plot(sol.t, sol.y[:2, :].T)    # first two components of solution
xlabel("t"), ylabel("angle")
title("Uncoupled pendulums");

You can see that the pendulums swing independently:

Source

from matplotlib import animation
fig, ax = subplots(1, 2)
rod_ = [] 
bob_ = []
for k in [0, 1]:
    rod_.append(ax[k].plot([], [], linewidth=3)[0])
    bob_.append(ax[k].scatter([0], [0], color='black', s=50))
    ax[k].set_aspect('equal')
    ax[k].set_xlim([-1.1, 1.1])
    ax[k].set_ylim([-1.1, 0.4])
    ax[k].set_axis_off()
time_ = ax[0].text(-0.95, 0.25, "", fontsize=11)

def snapshot(i):
    q = sol.y[:2, i]
    for k in [0, 1]:
        rod_[k].set_data([0, sin(q[k])], [0, -cos(q[k])])
        bob_[k].set_offsets([sin(q[k]), -cos(q[k])])
    time_.set_text(f"t = {sol.t[i]:.1f} s")

anim = animation.FuncAnimation(fig, snapshot, frames=401)
anim.save("figures/pendulums-weak.mp4", fps=24)

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 = 0.75    # changes the value in the du_dt closure
sol = solve_ivp(du_dt, tspan, u0, t_eval=linspace(0, 50, 1001))
plot(sol.t, sol.y[:2, :].T)
xlabel("t"), ylabel("angle")
title("Coupled pendulums");

The coupling makes the pendulums swap energy back and forth:

Source

anim = animation.FuncAnimation(fig, snapshot, frames=1001)
anim.save("figures/pendulums-strong.mp4", fps=24)

6.3.3Exercises¶

Exercise 6.3.1

✍ 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$

Exercise 6.3.2

✍ 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$

Exercise 6.3.3

⌨ Solve the following IVPs using Function 6.2.1 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]$

Exercise 6.3.5

⌨ 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]$ .

Exercise 6.3.6

⌨ 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