Systems and Higher-Order Equations

Every one-step method works on vectors as well as numbers. Any higher-order equation becomes a first-order system, as in Systems in Matrix Form, and then Euler, Heun or RK4 run with each $y$ replaced by a vector $\mathbf{u}$ and each $f$ by a vector function $\mathbf{F}$. The error orders stay the same.

Key idea

To solve $y'' = f(t, y, y')$, set $u_1 = y$, $u_2 = y'$: $$u_1' = u_2,\qquad u_2' = f(t, u_1, u_2).$$ Each component needs its own starting value, so the problem must give both $y(0)$ and $y'(0)$. Then step all components together, always using values from the same stage.

Why it works

The derivation of each method (Taylor expansion, or averaging slopes) never used that $y$ was a single number; it works component by component. What must not happen is mixing stages: in Euler's method both components of $\mathbf{u}_{n+1}$ are computed from the old $\mathbf{u}_n$. Updating $u_1$ first and then using the new $u_1$ in $u_2$ is a different method, with different behavior.

Worked example: Rewriting a second-order equation

Write $y'' + 0.5y' + 4y = 0$ as $\mathbf{u}' = A\mathbf{u}$.

$u_1 = y$, $u_2 = y'$, so $u_1' = u_2$.
$y'' = -4y - 0.5y'$, so $u_2' = -4u_1 - 0.5u_2$.
$A = \begin{bmatrix} 0 & 1 \\ -4 & -0.5 \end{bmatrix}$.

Answer: $\begin{bmatrix} 0 & 1 \\ -4 & -0.5 \end{bmatrix}$

Worked example: Euler on an oscillator

For $y'' = -y$, $y(0) = 1$, $y'(0) = 0$, take two Euler steps with $h = 0.1$.

System: $y' = v$, $v' = -y$, starting at $(1, 0)$.
Step 1: $y_1 = 1 + 0.1(0) = 1$, $v_1 = 0 + 0.1(-1) = -0.1$.
Step 2: $y_2 = 1 + 0.1(-0.1) = 0.99$, $v_2 = -0.1 + 0.1(-1) = -0.2$.
The exact solution keeps $y^{2} + v^{2} = 1$; Euler gives $1$, then $1.01$, then $1.0201$: it multiplies by $1 + h^{2}$ every step, spiraling slowly outward.

Answer: $(0.99, -0.2)$

Worked example: Heun does better

Take one Heun step for the same system from $(1, 0)$ with $h = 0.1$.

$\mathbf{k}_1 = \mathbf{F}(1, 0) = (0, -1)$. Predictor: $(1, 0) + 0.1(0, -1) = (1, -0.1)$.
$\mathbf{k}_2 = \mathbf{F}(1, -0.1) = (-0.1, -1)$.
$\mathbf{u}_1 = (1, 0) + 0.05\big((0, -1) + (-0.1, -1)\big) = (0.995, -0.1)$.
Exact: $(\cos 0.1, -\sin 0.1) \approx (0.995004, -0.099833)$.

Answer: $(0.995, -0.1)$

The system from the first example: trajectories spiral into the origin as the oscillation dies out. A good numerical method must follow these spirals; Euler with a large step spirals outward instead.

In code, $\mathbf{u}$ is simply an array, and $\mathbf{F}$ returns an array of the same length; the stepping loop does not change at all. Systems also show why stability matters: the components can decay or oscillate at very different rates, set by the eigenvalues of $A$, and a step size must suit all of them. That is the subject of Stability of Euler's Method.

Common mistake

Forgetting that each component needs its own starting value. A second-order equation needs both $y(0)$ and $y'(0)$; an $n$th-order equation needs $n$ values. Without $y'(0)$ the first step of $u_1' = u_2$ cannot even begin. Also keep the stages straight: Euler uses the old values for every component.

Sources and further reading

Tobin A. Driscoll and Richard J. Braun, Fundamentals of Numerical Computation (University of Delaware, free online edition): chapter "Initial-value problems for ODEs", section "IVP systems"

Autar Kaw, Holistic Numerical Methods (University of South Florida): Ordinary Differential Equations, "Runge–Kutta 2nd Order Method" (higher-order equations as systems)

Systems and Higher-Order Equations

Practice

Builds on

Leads to

Sources and further reading