Chapter 6

The OrdinaryDiffEq package offers solvers for IVPs. Let’s use it to define and solve an initial-value problem for $u'=\sin[(u+t)^2]$ over $t \in [0,4]$, such that $u(0)=-1$.

Because many practical problems come with parameters that are fixed within an instance but varied from one instance to another, the syntax for IVPs includes a input argument p that stays fixed throughout the solution. Here we don’t want to use that argument, but it must be in the definition for the solver to work.

f(u, p, t) = sin((t + u)^2)     # defines du/dt, must include p argument
u₀ = -1.0                       # initial value
tspan = (0.0, 4.0)               # t interval

With the data above we define an IVP problem object and then solve it. Here we tell the solver to use the Tsit5 method, which is a good first choice for most problems.

using OrdinaryDiffEq
ivp = ODEProblem(f, u₀, tspan)
sol = solve(ivp, Tsit5());

The resulting solution object can be shown using plot.

using Plots
plot(sol;
    label="solution", legend=:bottom,
    xlabel="t",  ylabel=L"u(t)",
    title=L"u'=\sin((t+u)^2)")

The solution also acts like any callable function that can be evaluated at different values of $t$.

@show sol(1.0);

sol(1.0) = -0.7903205813665345

Under the hood, the solution object holds some information about how the values and plot are produced:

[sol.t sol.u]

The solver initially finds approximate values of the solution (second column above) at some automatically chosen times (first column above). To compute the solution at other times, the object performs an interpolation on those values. This chapter is about how the discrete $t$ and $u$ values are computed. For now, just note how we can extract them from the solution object.

scatter!(sol.t, sol.u, label="discrete values")

The equation $u'=(u+t)^2$ gives us some trouble.

f(u, p, t) = (t + u)^2
ivp = ODEProblem(f, 1.0, (0.0, 1.0))
sol = solve(ivp, Tsit5());

┌ Warning: At t=0.7853839417697203, dt was forced below floating point epsilon 1.1102230246251565e-16, and step error estimate = 42.16290621054672. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase ~/.julia/packages/SciMLBase/XzPx0/src/integrator_interface.jl:623

The warning message we received can mean that there is a bug in the formulation of the problem. But if everything has been done correctly, it suggests that the solution may not exist past the indicated time. This is a possibility in nonlinear ODEs.

plot(sol, label="";
    xlabel=L"t",  yaxis=(:log10, L"u(t)"),
    title="Finite-time blowup")

Consider the ODEs $u'=u$ and $u'=-u$. In each case we compute $\partial f/\partial u = \pm 1$, so the condition number bound from Theorem 6.1.2 is $e^{b-a}$ in both problems. However, they behave quite differently. In the case of exponential growth, $u'=u$, the bound is the actual condition number.

t = range(0, 3, length=800)
u = @. exp(t) * 1
lower, upper = @. exp(t) * 0.7, @. exp(t) * 1.3
plot(t, u;
    l=:black, ribbon=(lower, upper),
    leg=:none,  xlabel=L"t",  ylabel=L"u(t)",
    title="Exponential divergence of solutions")

But with $u'=-u$, solutions actually get closer together with time.

u = @. exp(-t) * 1
lower, upper = @. exp(-t) * 0.7, @. exp(-t) * 1.3
plot(t, u;
    l=:black,  ribbon=(lower, upper),
    leg=:none,  xlabel=L"t",  ylabel=L"u(t)",
    title="Exponential convergence of solutions")

In this case the actual condition number is one, because the initial difference between solutions is the largest over all time. Hence the exponentially growing bound $e^{b-a}$ is a gross overestimate.

We consider the IVP $u'=\sin[(u+t)^2]$ over $0\le t \le 4$, with $u(0)=-1$.

using OrdinaryDiffEq
f(u, p, t) = sin((t + u)^2);
tspan = (0.0, 4.0);
u0 = -1.0;
ivp = ODEProblem(f, u0, tspan)

Here is the call to Function 6.2.2.

using Plots
t, u = FNC.euler(ivp, 20)
plot(t, u;
    m=2,  label="n=20", 
    xlabel=L"t",  ylabel=L"u(t)",
    title="Solution by Euler's method")

We could define a different interpolant to get a smoother picture above, but the derivation of Euler’s method assumed a piecewise linear interpolant. We can instead request more steps to make the interpolant look smoother.

t, u = FNC.euler(ivp, 50)
plot!(t, u, m=2, label="n=50")

Increasing $n$ changed the solution noticeably. Since we know that interpolants and finite differences become more accurate as $h\to 0$, we should anticipate the same behavior from Euler’s method. We don’t have an exact solution to compare to, so we will use a DifferentialEquations solver to construct an accurate reference solution.

u_exact = solve(ivp, Tsit5(), reltol=1e-14, abstol=1e-14)
plot!(u_exact, l=(2, :black), label="reference")

Now we can perform a convergence study.

n = [round(Int, 5 * 10^k) for k in 0:0.5:3]
err = []
for n in n
    t, u = FNC.euler(ivp, n)
    push!(err, norm(u_exact.(t) - u, Inf))
end
@pt :header=["n", "inf-norm error"] [n err]

The error is approximately cut by a factor of 10 for each increase in $n$ by the same factor. A log-log plot also confirms first-order convergence. Keep in mind that since $h=(b-a)/n$, it follows that $O(h)=O(n^{-1})$.

plot(n, err;
    m=:o, label="results",
    xaxis=(:log10, L"n"),  yaxis=(:log10, "inf-norm global error"),
    title="Convergence of Euler's method")

# Add line for perfect 1st order.
plot!(n, 0.5 * err[end] * (n / n[end]) .^ (-1), l=:dash, label=L"O(n^{-1})")

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, [α, β])

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.

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.

We solve the IVP $u'=\sin[(u+t)^2]$ over $0\le t \le 4$, with $u(0)=-1$.

using OrdinaryDiffEq
f(u, p, t) = sin((t + u)^2)
tspan = (0.0, 4.0)
u₀ = -1.0
ivp = ODEProblem(f, u₀, tspan)

We use a DifferentialEquations solver to construct an accurate approximation to the exact solution.

u_ref = solve(ivp, Tsit5(), reltol=1e-14, abstol=1e-14);

Now we perform a convergence study of our two Runge–Kutta implementations.

n = [round(Int, 2 * 10^k) for k in 0:0.5:3]
err = zeros(length(n), 2)
for (k, n) in enumerate(n)
    t, u = FNC.ie2(ivp, n)
    err[k, 1] = norm(u_ref.(t) - u, Inf)
    t, u = FNC.rk4(ivp, n)
    err[k, 2] = norm(u_ref.(t) - u, Inf)
end
@pt :header=["n", "IE2 error", "RK4 error"] [n err]

The amount of computational work at each time step is assumed to be proportional to the number of stages. Let’s compare on an apples-to-apples basis by using the number of $f$-evaluations on the horizontal axis.

using Plots
plot([2n 4n], err;
    m=3, label=["IE2" "RK4"], legend=:bottomleft,
    xaxis=(:log10, "f-evaluations"),  yaxis=(:log10, "inf-norm error"),
    title="Convergence of RK methods")

plot!(2n, 0.1 * err[end,1] * (n / n[end]) .^ (-2), l=:dash, label=L"O(n^{-2})")
plot!(4n, 0.1 * err[end,2] * (n / n[end]) .^ (-4), l=:dash, label=L"O(n^{-4})")

The fourth-order variant is more efficient in this problem over a wide range of accuracy.

Let’s run adaptive RK on $u'=e^{t-u\sin u}$.

using OrdinaryDiffEq, Plots
f(u, p, t) = exp(t - u * sin(u))
ivp = ODEProblem(f, 0, (0.0, 5.0))
t, u = FNC.rk23(ivp, 1e-5)
plot(t, u, m=2,
    xlabel=L"t",  ylabel=L"u(t)", 
    title="Adaptive IVP solution")

The solution makes a very abrupt change near $t=2.4$. The resulting time steps vary over three orders of magnitude.

Δt = diff(t)
plot(t[1:end-1], Δt;
    xaxis=(L"t", (0, 5)), yaxis=(:log10, "step size"),
    title="Adaptive step sizes")

If we had to run with a uniform step size to get this accuracy, it would be

println("minimum step size = $(minimum(Δt))")

minimum step size = 4.6096854609878335e-5

On the other hand, the average step size that was actually taken was

println("average step size = $(sum(Δt)/(length(t)-1))")

average step size = 0.03205128205128205

We took fewer steps by a factor of almost 1000! Even accounting for the extra stage per step and the occasional rejected step, the savings are clear.

In Demo 6.1.3 we saw an IVP that appears to blow up in a finite amount of time. Because the solution increases so rapidly as it approaches the blowup, adaptive stepping is required even to get close.

f(u, p, t) = (t + u)^2
ivp = ODEProblem(f, 1, (0.0, 1.0))
t, u = FNC.rk23(ivp, 1e-5);

┌ Warning: Stepsize too small near t=0.7854087204072808
└ @ FNCFunctions ~/Documents/GitHub/FNCFunctions.jl/src/chapter06.jl:102

In fact, the failure of the adaptivity gives a decent idea of when the singularity occurs.

plot(t, u;
    legend=:none,
    xlabel=L"t",  yaxis=(:log10, L"u(t)"), 
    title="Finite-time blowup")

tf = t[end]
vline!([tf], l=:dash)
annotate!(tf, 1e5, latexstring(@sprintf("t = %.6f ", tf)), :right)

We study the convergence of AB4 using the IVP $u'=\sin[(u+t)^2]$ over $0\le t \le 4$, with $u(0)=-1$. As usual, solve is called to give an accurate reference solution.

using OrdinaryDiffEq
ivp = ODEProblem((u, p, t) -> sin((t + u)^2), -1.0, (0.0, 4.0))
u_ref = solve(ivp, Tsit5(), reltol=1e-14, abstol=1e-14);

Now we perform a convergence study of the AB4 code.

n = @. [round(Int, 4 * 10^k) for k in 0:0.5:3]
err = []
for n in n
    t, u = FNC.ab4(ivp, n)
    push!(err, norm(u_ref.(t) - u, Inf))
end
@pt :header=["n", "inf-norm error"] [n err]

The method should converge as $O(h^4)$, so a log-log scale is appropriate for the errors.

using Plots
plot(n, err, m=3, 
    label="AB4",  legend=:bottomleft,
    xaxis=(:log10, L"n"),  yaxis=(:log10, "inf-norm error"),
    title="Convergence of AB4")

plot!(n, 0.1 * err[end] * (n / n[end]) .^ (-4), l=:dash, label=L"O(n^{-4})")

The following simple ODE uncovers a surprise.

ivp = ODEProblem((u, p, t) -> u^2 - u^3, 0.005, (0, 400.0))

We will solve the problem first with the implicit AM2 method using $n=200$ steps.

tI, uI = FNC.am2(ivp, 200)

plot(tI, uI;
    label="AM2", legend=:bottomright,
    xlabel=L"t",  ylabel=L"u(t)")

Now we repeat the process using the explicit AB4 method.

tE, uE = FNC.ab4(ivp, 200)
scatter!(tE, uE, m=3, label="AB4", ylim=[-4, 2])

Once the solution starts to take off, the AB4 result goes catastrophically wrong.

uE[105:111]

We hope that AB4 will converge in the limit $h\to 0$, so let’s try using more steps.

plt = scatter(tI, uI;
    m=3,  label="AM2, n=200",  legend=:bottomright,
    xlabel=L"t",  ylabel=L"u(t)")

for n in [1000, 1600]
    tE, uE = FNC.ab4(ivp, n)
    plot!(tE, uE, label="AM4, n=$n")
end
plt

So AB4, which is supposed to be more accurate than AM2, actually needs something like 8 times as many steps to get a reasonable-looking answer!

We’ll measure the error at the time $t=1$.

du_dt(u, t) = u
û = exp
a, b = 0.0, 1.0;
n = [5, 10, 20, 40, 60]
err = []
t, u = [], []
for n in n
    h = (b - a) / n
    t = [a + i * h for i in 0:n]
    u = [1; û(h); zeros(n - 1)]
    f_val = [du_dt(u[1], t[1]); zeros(n)]
    for i in 2:n
        f_val[i] = du_dt(u[i], t[i])
        u[i+1] = -4 * u[i] + 5 * u[i-1] + h * (4 * f_val[i] + 2 * f_val[i-1])
    end
    push!(err, abs(û(b) - u[end]))
end
@pt :header=["n", "h", "error"] [n (b - a) ./ n err]

The error starts out promisingly, but things explode from there. A graph of the last numerical attempt yields a clue.

using Plots
plot(t, abs.(u);
    m=3,  label="",
    xlabel=L"t",  yaxis=(:log10, L"|u(t)|"), 
    title="LIAF solution")

It’s clear that the solution is growing exponentially in time.