Vector-valued IVP

We can pack the two dependent variables \(y\) and \(z\) into a vector-valued function of time, \(\mathbf{u}(t)\), as \[\begin{align*} 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{align*}\] where \(u_1=y\), \(u_2=z\).

Numerical solutions

Generalizing Euler’s method to IVP systems is fairly straightforward, \[\mathbf{u}_{i+1} = \mathbf{u}_i + h\,\mathbf{f}(t_i,\mathbf{u}_i), \qquad i=0,\cdots,n-1.\]

"""
    euler(ivp,n)

Apply Euler's method to solve the given IVP using 'n' time steps.
Returns a vector of times and a vector of solution values.
"""
function euler(ivp,n)
        # Time discretization
        a,b = ivp.tspan
        h = (b-a)/n 
        t = [ a + i*h for i in 0:n ]

        # Initial condition and output setup.
        u = fill(float(ivp.u0),n+1)

        # The time stepping criterion
        for i in 1:n 
            u[i+1] = u[i] + h*ivp.f(u[i],ivp.p,t[i])
        end
        return t,u 
    end

euler

When there are two solutions, we can plot the solution in the phase plane.

# solution approaches a periodic form
plot(sol,vars=(1,2),title="Predator-prey in the phase plane",
    xlabel="y",ylabel="z")

Transformation of high-order systems

First-order ODE systems yields an approach for systems of higher-order as well.

With these definitions the system becomes \[\begin{align*} 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{align*}\] 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)\).

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.,50.);

γ,L,k = 0,0.5,0   # k = 0 => uncoupled
ivp = ODEProblem(couple,u₀,tspan,[γ,L,k])
sol = solve(ivp,Tsit5());

With coupling, we see the pendulums in phase as they are passing energy back and forth.

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

A second-order method

Consider a series expansion of the exact solution to \(u'=f(t,u)\), \[\hat{u}(t_{i+1}) = \hat{u}(t_i) + h \hat{u}'(t_i) + \frac{1}{2}h^2 \hat{u}''(t_i) + O(h^3).\]

Using this expression, we obtain \[\hat{u}(t_{i+1}) = \hat{u}(t_i) + h\left[f\bigl(t_i,\hat{u}(t_i)\bigr) + \frac{h}{2}f_t\bigl(t_i,\hat{u}(t_i)\bigr) + \frac{h}{2}f\bigl(t_i,\hat{u}(t_i)\bigr)\,f_u\bigl(t_i,\hat{u}(t_i)\bigr)\right] + O(h^3).\]

Implementation

Runge–Kutta methods are called multistage methods. In the first stage, the method takes an Euler half-step to time \(t_i+\frac{1}{2}h\):

\[\begin{align*} k_1 &= h f(t_i,u_i), \\ v &= u_i + \tfrac{1}{2}k_1. \end{align*}\]

Other Runge-Kutta methods

Using the Taylor series derivation, one can derive other methods of even higher orders of accuracy. Doing so requires extra stages, and the algebra grows rapidly in complexity. We won’t go into details regarding this derivation.

This recipe is completely determined by the number of stages \(s\) and the constants \(a_{ij}\), \(b_j\), and \(c_i\). Often an RK method is presented as just a table of these numbers, as in

\[\begin{array}{r|ccccc} 0 & & & & & \\ c_1 & a_{11} & & &\\ c_2 & a_{21} & a_{22} & & &\\ \vdots & \vdots & & \ddots & &\\ c_{s-1} & a_{s-1,1} & \cdots & & a_{s-1,s-1}&\\[1mm] \hline \rule{0pt}{2.25ex} & b_1 & b_2 & \cdots & b_{s-1} & b_s \end{array}\]

Here are two more two-stage, second-order methods, modified Euler and Heun’s method, respectively:

\[\begin{array}{r|cc} \rule{0pt}{2.75ex}0 & & \\ \rule{0pt}{2.75ex}1 & 1 &\\[1mm] \hline \rule{0pt}{2.75ex}& \frac{1}{2} & \frac{1}{2} \end{array}\quad\quad \begin{array}{r|cc} \rule{0pt}{2.75ex} 0 & & \\ \rule{0pt}{2.75ex} \frac{2}{3} & \frac{2}{3} &\\[1mm] \hline \rule{0pt}{2.75ex} & \frac{1}{4} & \frac{3}{4} \end{array}\]

Fourth-order Runge-Kutta

The most commonly used RK method, and perhaps the most popular IVP method of all, is the fourth-order one given by

\[\begin{array}{r|cccc} \rule{0pt}{2.75ex}0 & & & & \\ \rule{0pt}{2.75ex}\frac{1}{2} & \frac{1}{2} & & &\\ \rule{0pt}{2.75ex}\frac{1}{2} & 0 & \frac{1}{2} & &\\ \rule{0pt}{2.75ex}1 & 0 & 0 & 1\\[1mm] \hline \rule{0pt}{2.75ex}& \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array}\]