14.1. Initial Value Problems

An initial value problem (IVP) can be written as:

\[\begin{split} \begin{align*} y'(t) &= f(t,y(t)), \quad t\ge t_0 \\ y(t_0) &= y_0 \end{align*} \end{split}\]

for some right-hand side function \(f(t,y)\) and initial condition \(y_0\) for \(t=t_0\). In these problems, the independent variable \(t\) often represents time.

Sometimes basic calculus can be used to solve these equations. For example, some linear and separable equations have closed-form solutions. But in general we have to resort to numerical methods to find approximations to the solution \(y(t)\).

using PyPlot, PyCall

14.1.1. Euler’s method

The most basic approach for numerically solving an IVP is Euler’s method, illustrated below.

(from Wikipedia, https://en.wikipedia.org/wiki/Euler_method).

We first choose a step size, or time step, \(h\), and define the discrete times \(t_n = t_0 + nh\), \(n = 0,1,\ldots\). Euler’s method advances the solution \(y_n\) at time \(t_n\) to time \(t_{n+1}\) by a linear approximation using the derivative at time \(t_n\):

\[ y_{n+1} = y_n + hf(t_n,y_n) \]

With some assumptions on the IVP and the step size, these values can be shown to be good approximation to the true solution, that is, \(y_n \approx y(t_n\)).

While described for scalar-valued functions \(y(t)\), the method works equally well for systems of equations, that is, when \(y(t)\) and \(f(t,y(t))\) are vector-valued.

14.1.2. Euler’s method, implementation

We implement a general version of Euler’s method which takes the right-hand side function \(f(t,y(t))\) as an argument. In addition, it takes the initial condition \(y_0\), the step size \(h\), the number of steps \(N\), and the starting time \(t_0\) which has a default value of \(0\). Note that the implementations supports solutions \(y\) that are vector-valued, and the output y is a 2D array of the approximate solutions at each time step.

function euler(f, y0, h, N, t0=0.0)
    t = t0 .+ h*(0:N)
    y = zeros(N+1, length(y0))
    
    y[1,:] .= y0
    for n = 1:N
        y[n+1,:] = y[n,:] + h * f(t[n], y[n,:])
    end
    
    return t,y
end

euler (generic function with 2 methods)

We demonstrate the method on a model problem with \(f(t,y) = -y + \sin t\) and \(y(0) = 1\). We solve using Euler’s method with time step \(h=0.2\) and \(N=20\) steps.

f(t,y) = -y .+ sin(t)
t,y = euler(f, 1, 0.2, 20)
plot(t, y, "-o");

We can compare the results to the exact solution:

\[ y_\mathrm{exact} = e^{-t} + (\sin t - \cos t + e^{-t}) / 2 \]

In the code below, we solve three times using various time steps \(h\). In the plot, we see that the results are generally better for smaller \(h\) (as expected), and the accuracy seems to roughly scale linearly with \(h\). This can be shown to be true in general: the method is convergent (the approximations approach the exact solution as \(h \rightarrow 0\)), and the method is first-order accurate (the error scales with the first power of \(h\)).

yexact(t) = exp(-t) + (sin(t) - cos(t) + exp(-t)) / 2
tt = 0:0.01:2
plot(tt, yexact.(tt))
for h = [0.5, 0.2, 0.1]
    t,y = euler(f, 1, h, round(Int, 2/h))
    plot(t, y, "-o")
end
legend(("Exact", "h=0.5", "h=0.2", "h=0.1"));

14.1.3. The Runge-Kutta method

While Euler’s method does give accurate results for small enough time steps \(h\), its first-order convergence is often considered too slow and much better methods have been developed. One of the most popular methods is the following 4th order accurate Runge-Kutta method (RK4):

\[\begin{split} \begin{align*} k_1 &= h f(t_n, y_n) \\ k_2 &= h f(t_n + h/2, y_n + k_1/2) \\ k_3 &= h f(t_n + h/2, y_n + k_2/2) \\ k_4 &= h f(t_n + h, y_n + k_3) \\ y_{n+1} &= y_n + (k_1 + 2k_2 + 2k_3 + k_4) / 6 \end{align*} \end{split}\]

We implement it in a very similar way to before. Note that each step now requires computing the four so-called stages or stage derivatives \(k_1,k_2,k_3,k_4\).

function rk4(f, y0, h, N, t0=0)
    t = t0 .+ h*(0:N)
    y = zeros(N+1, length(y0))
    
    y[1,:] .= y0
    for n = 1:N
        k1 = h * f(t[n], y[n,:])
        k2 = h * f(t[n] + h/2, y[n,:] + k1/2)
        k3 = h * f(t[n] + h/2, y[n,:] + k2/2)
        k4 = h * f(t[n] + h, y[n,:] + k3)
        y[n+1,:] = y[n,:] + (k1 + 2k2 + 2k3 + k4) / 6
    end
    
    return t,y
end

rk4 (generic function with 2 methods)

Run the same test problem as before using RK4:

yexact(t) = exp(-t) + (sin(t) - cos(t) + exp(-t)) / 2
tt = 0:0.01:2
plot(tt, yexact.(tt))
for h = [0.5, 0.2, 0.1]
    t,y = rk4(f, 1, h, round(Int, 2/h))
    plot(t, y, "-o")
end
legend(("Exact", "h=0.5", "h=0.2", "h=0.1"));

We can see that the results are fundamentally more accurate. In fact, even with largest time step \(h=0.5\), the point values are essentially right on top of the exact solution. We can compute the error for the smallest time step \(h=0.1\), and note that it has 6 accurate digits:

errors = @. y - yexact(t)
max_error = maximum(abs.(errors))

5.015261516083669e-7