Chapter 6

The DifferentialEquations 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$.

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.

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

The resulting solution object can be shown using plot.

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);

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());

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

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.

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

pretty_table((n=n, err=err), header=["n", "Inf-norm error"])

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.05 * (n / n[1]) .^ (-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.

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.

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];

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.

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.

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);

γ, 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$.

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_IE2, err_RK4 = [], []
for n in n
    t, u = FNC.ie2(ivp, n)
    push!(err_IE2, maximum(@.abs(u_ref(t) - u)))
    t, u = FNC.rk4(ivp, n)
    push!(err_RK4, maximum(@.abs(u_ref(t) - u)))
end

pretty_table([n err_IE2 err_RK4], header=["n", "IE2 error", "RK4 error"])

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.

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

plot!(2n, 1e-5 * (n / n[end]) .^ (-2), l=:dash, label=L"O(n^{-2})")
plot!(4n, 1e-10 * (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}$.

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, title="Adaptive step sizes",
    xaxis=(L"t", (0, 5)), yaxis=(:log10, "step size"))

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

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

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

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

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);

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.

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

pretty_table([n err], header=["n", "inf-norm error"])

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

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

plot!(n, (n / n[1]) .^ (-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",
    xlabel=L"t", ylabel=L"u(t)", leg=:bottomright)

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, label="AM2, n=200", m=3,
    xlabel=L"t", ylabel=L"u(t)", leg=:bottomright)

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

pretty_table([n (b - a) ./ n err], header=["n", "h", "error"])

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

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.