Chapter 5 — Initial-Value Problems for Ordinary
Differential Equations
Math 128A, Fall 2026 · Per-Olof Persson
Most differential equations that describe something real cannot be
solved in closed form, so we solve them numerically: start from the
initial condition and step forward. This chapter is about how to take
that step, and about what can go wrong. Two themes recur. Order
says how fast the error shrinks as the step size does, and higher is
better. Stability says whether the method’s own errors stay
bounded, and it can override order entirely — a high-order method used
outside its stability region produces garbage, however small the
step.
5.1 The
Elementary Theory of Initial-Value Problems
Where these problems come
from
Example 5.1
The pendulum.
A mass \(m\) hangs on a rigid rod of
length \(\ell\), displaced by an angle
\(\theta\) from the vertical.
Figure 5.1: The pendulum.
Gravity’s tangential component drives the motion; the radial component
is taken up by the rod.
Newton’s second law is \(\vec F = m\vec
a\). Along the arc, the tangential component of gravity is \[
F_t = -mg\sin\theta,
\] the minus sign because it pulls back towards the vertical, and
the tangential acceleration is \(a =
\ell\theta''\). So \[
-mg\sin\theta = m\ell\theta''
\qquad\Longrightarrow\qquad
\theta'' + \frac{g}{\ell}\sin\theta = 0,
\] a nonlinear second-order ordinary differential
equation.
For small \(\theta\) we may use
\(\sin\theta \approx \theta\), giving
the linear equation \[
\theta'' + \frac{g}{\ell}\theta = 0 .
\] Either way the equation alone does not determine the motion:
we also need the initial state, \[
\theta(t_0) = \theta_0, \qquad \theta'(t_0) = \omega_0,
\] which is what makes this an initial-value
problem. With \(\theta(0)=\theta_0\) and \(\theta'(0)=0\) the linearised problem
has the solution \[
\theta(t) = \theta_0\cos\!\left(\sqrt{\tfrac{g}{\ell}}\,t\right).
\]
Three degrees of smoothness
The standard first-order initial-value problem is \[
y'(t) = f(t, y(t)), \qquad a \le t \le b, \qquad y(a) = \alpha .
\]
Everything in the chapter depends on how smooth \(f\) is in its second argument, and three
conditions need distinguishing:
Lipschitz continuity — \(|f(x) - f(y)| \le \lambda |x-y|\) for a
fixed constant \(\lambda\).
Differentiability — \(\dfrac{|f(x)-f(y)|}{|x-y|} \to \mu\) as
\(|x-y| \to 0\).
Each is strictly stronger than the one before, and Lipschitz
continuity — the middle one — is exactly what the existence and
uniqueness theory needs.
Lipschitz Condition and Convexity
Definition 5.1
A function \(f(t,y)\) is said to
satisfy a Lipschitz condition in the variable \(y\) on a set \(D\subset \mathbb{R}^2\) if a constant \(L>0\) exists with \[
\begin{aligned}
|f(t,y_1)-f(t,y_2)|\le L |y_1-y_2|,
\end{aligned}
\] whenever \((t,y_1),(t,y_2)\in
D\). The constant \(L\) is
called a Lipschitz constant for \(f\).
Definition 5.2
A set \(D\subset \mathbb{R}^2\) is
said to be convex if whenever \((t_1,y_1)\) and \((t_2,y_2)\) belong to \(D\) and \(\lambda\) is in \([0,1]\), the point \(((1-\lambda)t_1+\lambda t_2,(1-\lambda)y_1+\lambda
y_2)\) also belongs to \(D\).
Example 5.2
The two conditions really are different.
\(y(x) = |x|\) is continuous, and
Lipschitz with \(\lambda = 1\), but it
is not differentiable at \(x
= 0\). So Lipschitz does not imply differentiable.
\(y(x) = \sqrt{x}\) is continuous on
\([0,\infty)\) but not
Lipschitz there: the difference quotient at the origin, \[
\frac{f(\varepsilon) - f(0)}{\varepsilon - 0} =
\frac{\sqrt\varepsilon}{\varepsilon}
= \frac{1}{\sqrt{\varepsilon}} \longrightarrow \infty
\qquad\text{as } \varepsilon \to 0,
\] is unbounded, so no single \(\lambda\) can work. It is not
differentiable at \(0\) either.
Figure 5.2: Left: \(|x|\) is Lipschitz but not differentiable
at the origin. Right: \(\sqrt{x}\) is
continuous but not Lipschitz on \([0,\infty)\).
Existence and uniqueness
Existence and Uniqueness
Theorem 5.1
Suppose \(f(t,y)\) is defined on a
convex set \(D\subset \mathbb{R}^2\).
If a constant \(L>0\) exists with
\[
\begin{aligned}
\left|\frac{\partial f}{\partial y}(t,y)\right|\le L,\quad
\text{for all }(t,y)\in D,
\end{aligned}
\] then \(f\) satisfies a
Lipschitz condition on \(D\) in the
variable \(y\) with Lipschitz constant
\(L\).
Theorem 5.2
Suppose that \(D=\{(t,y)\ |\ a\le t \le b,
-\infty < y < \infty \}\) and that \(f(t,y)\) is continuous on \(D\). If \(f\) satisfies a Lipschitz condition on
\(D\) in the variable \(y\), then the initial-value problem \[
\begin{aligned}
y'(t)=f(t,y),\quad a\le t\le b,\quad y(a)=\alpha,
\end{aligned}
\] has a unique solution \(y(t)\) for \(a\le
t\le b\).
In practice the Lipschitz condition is checked through the partial
derivative: if \(\partial f/\partial
y\) is bounded by \(L\) on the
region of interest, then \(f\) is
Lipschitz in \(y\) with constant \(L\).
Example 5.3
A well-behaved problem.
For \(y' = y\) we have \(f(t,y) = y\), so \(\partial f/\partial y = 1\) and \(f\) is Lipschitz with \(L = 1\). The theorem therefore promises a
unique solution, and indeed \(y(t) =
y_0e^{t}\).
Example 5.4
Existence can fail in finite time.
Consider \(y' = y^2\) with \(y(0) = y_0 > 0\). Here \[
\frac{\partial f}{\partial y} = 2y,
\] which is unbounded on \(D = \{(t,y)
\mid a \le t \le b,\ -\infty < y < \infty\}\), so \(f\) is not Lipschitz on
all of \(D\). Separating variables
gives \[
y(t) = \frac{1}{\frac{1}{y_0} - t},
\] which blows up as \(t \to
1/y_0\).
Figure 5.3: The solution of
\(y' = y^2\) reaches infinity at
\(t = 1/y_0\): no solution exists on
the whole of \(t \ge 0\).
So the solution genuinely does not exist for all \(t \ge 0\) — the theorem’s hypothesis was
not a technicality.
Example 5.5
Uniqueness can fail too.
Consider \(y' = 2\sqrt{y}\) for
\(t \ge 0\) with \(y(0) = 0\). As we saw, \(\sqrt{\cdot}\) is not Lipschitz near zero,
and the consequence is that the problem has infinitely
many solutions:
\(y(t) = 0\) for all \(t\);
\(y(t) = t^2\);
for any \(a \ge 0\), \[
y(t) = \begin{cases} 0, & t \le a, \\ (t-a)^2, & t > a.
\end{cases}
\]
Figure 5.4: Three of the
infinitely many solutions of \(y'=2\sqrt
y\), \(y(0)=0\). The solution
can stay at zero for any length of time before leaving it.
Each of these satisfies the equation everywhere, including at the
corner. A numerical method given this problem will return one
of them, and which one depends on the arithmetic — which is why the
hypothesis matters.
Example 5.6
A textbook check.
Show that \[
y' = \frac{2}{t}y + t^2e^t,
\qquad 1 \le t \le 2, \qquad y(1) = 0
\] has a unique solution.
Here \(f(t,y) = \dfrac{2}{t}y +
t^2e^t\), so \[
\frac{\partial f}{\partial y} = \frac{2}{t} \le 2
\qquad\text{for } 1 \le t \le 2 .
\] So \(f\) is Lipschitz in
\(y\) with \(L = 2\) on \(D =
\{(t,y) \mid 1 \le t \le 2,\ -\infty < y < \infty\}\), and
\(f\) is continuous on \(D\). The theorem applies: the solution
exists and is unique.
Well-Posedness
Definition 5.3
The initial-value problem \[
\begin{aligned}
\frac{dy}{dt} = f(t,y),\quad a\le t\le b,\quad y(a)=\alpha,
\end{aligned}
\] is said to be a well-posed problem if:
A unique solution, \(y(t)\), to the
problem exists, and
There exist constants \(\varepsilon_0>0\) and \(k>0\) such that for any \(\varepsilon\), with \(\varepsilon_0>\varepsilon>0\),
whenever \(\delta(t)\) is continuous
with \(|\delta(t)|<\varepsilon\) for
all \(t\) in \([a,b]\), and when \(|\delta_0|<\varepsilon\), the
initial-value problem \[
\begin{aligned}
\frac{dz}{dt}=f(t,z)+\delta(t),\quad a\le t \le b,\quad
z(a)=\alpha+\delta_0,
\end{aligned}
\] has a unique solution \(z(t)\) that satisfies \[
\begin{aligned}
|z(t)-y(t)|<k\varepsilon\quad\text{for all }t\text{ in }[a,b].
\end{aligned}
\]
Well-Posedness
Theorem 5.3
Suppose \(D=\{(t,y)\ |\ a\le t \le b\text{
and }-\infty<y<\infty\}\). If \(f\) is continuous and satisfies a Lipschitz
condition in the variable \(y\) on the
set \(D\), then the initial-value
problem \[
\begin{aligned}
\frac{dy}{dt} = f(t,y),\quad a\le t\le b,\quad y(a)=\alpha
\end{aligned}
\] is well-posed.
5.2 Euler’s Method
Euler’s Method
Suppose a well-posed initial-value problem is given: \[
\begin{aligned}
\frac{dy}{dt} = f(t,y),\quad a\le t \le b,\quad y(a)=\alpha
\end{aligned}
\] Distribute mesh points equally throughout \([a,b]\): \[
\begin{aligned}
t_i=a+ih,\quad\text{for each }i=0,1,2,\ldots,N.
\end{aligned}
\] The step size\(h=(b-a)/N =
t_{i+1}-t_i\).
Euler’s method is the simplest possible step: follow the tangent. Put
\[
h = \frac{b-a}{N},
\qquad t_i = a + ih,
\] and let \(w_i\) denote the
approximation to \(y(t_i)\). Since
\(y'(t_i) = f(t_i,
y(t_i))\) is the slope of the solution through that point,
stepping along it gives \[
w_0 = \alpha, \qquad w_{i+1} = w_i + h\,f(t_i, w_i).
\]
Figure 5.5: Euler’s method.
At each node the tangent slope \(f(t_i,w_i)\) is followed for one step of
\(h\). The polyline drifts away from
the solution as the steps accumulate.
Euler’s Method
Use Taylor’s Theorem for \(y(t)\):
\[
\begin{aligned}
y(t_{i+1}) = y(t_i) + (t_{i+1}-t_i)y'(t_i)
+\frac{(t_{i+1}-t_i)^2}{2}y''(\xi_i)
\end{aligned}
\] for \(\xi_i\in(t_i,t_{i+1})\). Since \(h=t_{i+1}-t_i\) and \(y'(t_i)=f(t_i,y(t_i))\), \[
\begin{aligned}
y(t_{i+1}) = y(t_i) + h f(t_i,y(t_i)) + \frac{h^2}{2}y''(\xi_i).
\end{aligned}
\] Neglecting the remainder term gives Euler’s method for \(w_i\approx y(t_i)\): \[
\begin{aligned}
w_0&=\alpha \\
w_{i+1} &= w_i + h f(t_i,w_i),\qquad i=0,1,\ldots,N-1
\end{aligned}
\] The well-posedness implies that \[
\begin{aligned}
f(t_i,w_i) \approx y'(t_i) = f(t_i,y(t_i))
\end{aligned}
\]
The error bound
Error Bound
Theorem 5.4
Suppose \(f\) is continuous and
satisfies a Lipschitz condition with constant \(L\) on \[
\begin{aligned}
D = \{(t,y)\ |\ a\le t \le b, -\infty<y<\infty \}
\end{aligned}
\] and that a constant \(M\)
exists with \[
\begin{aligned}
|y''(t)|\le M,\quad \text{for all }t\in[a,b].
\end{aligned}
\] Let \(y(t)\) denote the
unique solution to the initial-value problem \[
\begin{aligned}
y'=f(t,y),\quad a\le t \le b,\quad y(a)=\alpha,
\end{aligned}
\] and \(w_0,w_1,\ldots,w_n\) as
in Euler’s method. Then \[
\begin{aligned}
|y(t_i)-w_i|\le \frac{hM}{2L}\left[ e^{L(t_i-a)}-1\right].
\end{aligned}
\]
Proof of the Euler error bound.
Taylor expanding the solution about \(t_i\), \[
y(t_{i+1}) = y(t_i) + hy'(t_i) + \frac{h^2}{2}y''(\xi_i),
\] or, writing \(y_i = y(t_i)\)
and \(y_i' = y'(t_i)\), \[
y_{i+1} = y_i + hy_i' + \frac{h^2}{2}y''(\xi_i).
\]
Subtracting Euler’s method \(w_{i+1} = w_i
+ hf(t_i,w_i)\), \[
\begin{aligned}
y_{i+1} - w_{i+1}
&= y_i - w_i + h\bigl(y_i' - f(t_i,w_i)\bigr) +
\frac{h^2}{2}y''(\xi_i) \\
&= (y_i - w_i) + h\bigl(f(t_i,y_i) - f(t_i,w_i)\bigr) +
\frac{h^2}{2}y''(\xi_i),
\end{aligned}
\] using \(y_i' = f(t_i,
y_i)\). Taking absolute values and applying the Lipschitz
condition to the middle term and the bound \(|y''| \le M\) to the last, \[
|y_{i+1} - w_{i+1}|
\le |y_i - w_i| + hL|y_i - w_i| + \frac{h^2}{2}M
= (1 + hL)|y_i - w_i| + \frac{h^2}{2}M .
\]
This is a recurrence of the form \(a_{i+1}
= \gamma a_i + \beta\), for which an easy induction gives \[
a_i = \gamma^i a_0 + \beta\,\frac{\gamma^i - 1}{\gamma - 1}
\qquad (\gamma \ne 1).
\]
The bound is proportional to \(h\),
so in particular \(|y_i - w_i| \to 0\)
as \(h \to 0\): Euler’s method
converges, at first order. The exponential factor is a warning, though —
the bound grows with the length of the interval, and rapidly if \(L\) is large.
Local Truncation Error
Definition 5.4
The difference method \[
\begin{aligned}
w_0 &= \alpha \\
w_{i+1} &= w_i + h\phi(t_i,w_i)
\end{aligned}
\] has local truncation error\[
\begin{aligned}
\tau_{i+1}(h) =
\frac{y_{i+1}-(y_i+h\phi(t_i,y_i))}{h} =
\frac{y_{i+1}-y_i}{h} - \phi(t_i,y_i),
\end{aligned}
\] for each \(i=0,1,\ldots,N-1\).
5.3 Higher-Order Taylor
Methods
Higher-Order Taylor Methods
Consider initial-value problem \[
\begin{aligned}
y'=f(t,y),\quad a\le t \le b,\quad y(a)=\alpha.
\end{aligned}
\] Expand \(y(t)\) in \(n\)th Taylor polynomial about \(t_i\), evaluated at \(t_{i+1}\): \[
\begin{aligned}
y(t_{i+1}) &= y(t_i)+hy'(t_i)+\frac{h^2}{2}y''(t_i)
+ \cdots \\
&+\frac{h^n}{n!}y^{(n)}(t_i) + \frac{h^{n+1}}{(n+1)!}
y^{(n+1)}(\xi_i) \\
&= y(t_i) + hf(t_i,y(t_i)) + \frac{h^2}{2}f'(t_i,y(t_i)) +
\cdots \\
&+ \frac{h^n}{n!}f^{(n-1)}(t_i,y(t_i))+
\frac{h^{n+1}}{(n+1)!}f^{(n)}(\xi_i,y(\xi_i))
\end{aligned}
\] for some \(\xi_i\in
(t_i,t_{i+1})\). Delete remainder term to obtain the Taylor
method of order \(n\).
Euler’s method keeps one term of the Taylor expansion. Keeping more
gives higher order — the Taylor method of order \(n\): \[
w_{i+1} = w_i + hf(t_i,w_i) + \frac{h^2}{2}f'(t_i,w_i)
+ \cdots + \frac{h^n}{n!}f^{(n-1)}(t_i,w_i).
\]
The catch is that \(f'\) here
means the total derivative along the solution, so the chain
rule brings in partial derivatives of \(f\): \[
f'(t, y(t)) = \frac{\partial f}{\partial t}
+ \frac{\partial f}{\partial
y}\cdot\underbrace{\frac{dy}{dt}}_{=\,y'\,=\,f}
= f_t + f_y f .
\]
This problem has the exact solution \(y(t)
= y_0e^{t(1-8t)}\), which makes it convenient for checking a
method. The order-2 method is \[
w_0 = \alpha, \qquad
w_{i+1} = w_i + hf(t_i,w_i) + \frac{h^2}{2}f'(t_i,w_i),
\] so we need \(f'\). With
\(f(t,y(t)) = y(t)(1-16t)\), \[
\begin{aligned}
f'(t,y(t)) &= \frac{\partial f}{\partial t} + \frac{\partial
f}{\partial y}\cdot y' \\
&= -16y(t) + (1-16t)\,y(t)(1-16t) \\
&= -16y(t) + (1-16t)^2 y(t),
\end{aligned}
\] where the second term used \(y'
= f = y(1-16t)\). Substituting, \[
w_{i+1} = w_i + h\,f(t_i,w_i)
+ \frac{h^2}{2}\Bigl[-16w_i + (1-16t_i)^2 w_i\Bigr].
\]
Each extra order costs another derivative of \(f\), worked out by hand. That is the price
the next section is designed to avoid.
Higher-Order Taylor Methods
Theorem 5.5
If Taylor’s method of order \(n\) is
used to approximate the solution to \[
\begin{aligned}
y'(t) = f(t,y(t)),\quad a\le t\le b,\quad y(a)=\alpha,
\end{aligned}
\] with step size \(h\) and if
\(y\in C^{n+1}[a,b]\), then the local
truncation error is \(O(h^n)\).
5.4 Runge-Kutta Methods
Taylor’s Theorem in Two Variables
Theorem 5.6
Suppose \(f(t,y)\) and partial
derivatives up to order \(n+1\)
continuous on \(D=\{(t,y)\ |\ a\le t \le b,
c\le y \le d\}\), let \((t_0,y_0)\in
D\). For \((t,y)\in D\), there
is \(\xi\in[t,t_0]\) and \(\mu\in[y,y_0]\) with \[
\begin{aligned}
f(t,y) &= P_n(t,y)+ R_n(t,y) \\
P_n(t,y) &=f(t_0,y_0) + \left[(t-t_0)\frac{\partial f}{\partial t}
(t_0,y_0) + (y-y_0)\frac{\partial f}{\partial y}(t_0,y_0)\right] \\
&+\left[\frac{(t-t_0)^2}{2} \frac{\partial^2 f}{\partial
t^2}(t_0,y_0)
+(t-t_0)(y-y_0)\frac{\partial^2 f}{\partial t\partial y}
(t_0,y_0)\right. \\
&+\left.\frac{(y-y_0)^2}{2} \frac{\partial^2 f}{\partial y^2}
(t_0,y_0)\right] +
\cdots \\
&+ \left[ \frac{1}{n!} \sum_{j=0}^n \binom{n}{j} (t-t_0)^{n-j}
(y-y_0)^j \frac{\partial^n f}{\partial t^{n-j}\partial
y^j}(t_0,y_0)\right]
\end{aligned}
\]
Runge-Kutta methods get the order of a Taylor method without the
derivatives, by evaluating \(f\) at
cleverly chosen points instead.
Start from the order-2 Taylor increment, expanded with the chain
rule: \[
T^{(2)}(t,y) = f(t,y) + \frac{h}{2}f'(t,y)
= f(t,y) + \frac{h}{2}f_t(t,y) + \frac{h}{2}f_y(t,y)f(t,y).
\]
Can this be matched by a single evaluation of \(f\) at a shifted point, \[
a_1 f(t + \alpha_1,\, y + \beta_1)?
\]
Taylor’s Theorem in Two Variables
Theorem 5.7
(cont’d) \[
\begin{aligned}
R_n(t,y) &= \frac{1}{(n+1)!} \sum_{j=0}^{n+1}
\binom{n+1}{j} (t-t_0)^{n+1-j} (y-y_0)^j \cdot \\
&\cdot\frac{\partial^{n+1} f}{\partial t^{n+1-j}\partial
y^j}(\xi,\mu)
\end{aligned}
\]\(P_n(t,y)\) is the \(n\)th Taylor polynomial in two
variables.
Expanding that in a two-variable Taylor series, \[
a_1 f(t,y) + a_1\alpha_1 f_t(t,y) + a_1\beta_1 f_y(t,y) + a_1 R_1,
\]
and matching coefficients against \(T^{(2)}\) term by term: \[
a_1 = 1,
\qquad
\alpha_1 = \frac{h}{2},
\qquad
\beta_1 = \frac{h}{2}f(t,y).
\]
So \[
T^{(2)}(t,y) = f\!\left(t + \frac{h}{2},\, y + \frac{h}{2}f(t,y)\right)
+ O(h^2),
\] which gives the midpoint method: \[
w_0 = \alpha, \qquad
w_{i+1} = w_i + h\,f\!\left(t_i + \frac{h}{2},\, w_i +
\frac{h}{2}f(t_i,w_i)\right),
\] usually written in stages as \[
k_1 = h\,f(t_i, w_i),
\qquad
k_2 = h\,f\!\left(t_i + \frac{h}{2},\, w_i + \frac{k_1}{2}\right),
\qquad
w_{i+1} = w_i + k_2 .
\]
Figure 5.7: The midpoint
method. Euler’s step (blue) uses the slope at the left endpoint and
overshoots; the slope at the midpoint of that step (green) is a better
representative, and the step taken with it (magenta) lands much closer
to the solution.
Runge-Kutta Methods
Obtain high-order accuracy of Taylor methods without knowledges of
derivatives of \(f\)
Determine \(a_1,\alpha_1,\beta_1\)
such that \[
\begin{aligned}
a_1 f(t+\alpha_1,y+\beta_1) \approx
f(t,y) + \frac{h}{2} f'(t,y) = T^{(2)}(t,y).
\end{aligned}
\] with \(O(h^2)\) error.
Since \[
\begin{aligned}
f'(t,y) = \frac{df}{dt}(t,y) = \frac{\partial f}{\partial t}(t,y)
+\frac{\partial f}{\partial y}(t,y)\cdot y'(t)
\end{aligned}
\] and \(y'(t)=f(t,y)\), we
have \[
\begin{aligned}
T^{(2)}(t,y)=f(t,y)+\frac{h}{2}\frac{\partial f}{\partial t}(t,y)
+ \frac{h}{2}\frac{\partial f}{\partial y}(t,y)\cdot f(t,y)
\end{aligned}
\]
Averaging the two end slopes instead of using the midpoint gives
another second-order method, the explicit trapezoidal
(or modified Euler) method: \[
k_1 = h\,f(t_i, w_i),
\qquad
k_2 = h\,f(t_i + h,\, w_i + k_1),
\qquad
w_{i+1} = w_i + \frac{k_1 + k_2}{2}.
\]
Figure 5.8: The explicit
trapezoidal method. Euler’s slope predicts a point at \(t_1\); the slope there is averaged with the
original, and the step is taken with the average.
Runge-Kutta of order four is the one to remember. It costs four
evaluations of \(f\) per step and is
the default choice for smooth, non-stiff problems.
5.9 Higher-Order Equations
and Systems
High-Order Systems of Initial-Value Problems
An \(m\)th-order system of
first-order initial-value problems has the form \[
\begin{aligned}
\frac{du_1}{dt}(t) &= f_1(t,u_1,u_2,\ldots,u_m), \\
\frac{du_2}{dt}(t) &= f_2(t,u_1,u_2,\ldots,u_m), \\
& \vdots \\
\frac{du_m}{dt}(t) &= f_m(t,u_1,u_2,\ldots,u_m),
\end{aligned}
\] for \(a\le t \le b\), with
the initial conditions \[
\begin{aligned}
u_1(a)=\alpha_1, u_2(a)=\alpha_2,\ldots,u_m(a)=\alpha_m.
\end{aligned}
\]
Everything so far was written for a single scalar equation, and
nothing needs to change for a system: read \(y\), \(f\)
and \(\alpha\) as vectors. With \[
\vec y = \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_m \end{pmatrix},
\qquad
\vec f = \begin{pmatrix} f_1 \\ f_2 \\ \vdots \\ f_m \end{pmatrix},
\qquad
\vec\alpha = \begin{pmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_m
\end{pmatrix},
\] the problem is \(\vec y\,' =
\vec f(t, \vec y)\) with \(\vec y(a) =
\vec\alpha\), and Euler’s method reads \[
\vec w_{i+1} = \vec w_i + h\,\vec f(t_i, \vec w_i),
\qquad
\vec w_i = \begin{pmatrix} w_{1,i} \\ w_{2,i} \end{pmatrix}
\] for a two-component system — componentwise, the same
formula.
This is the reason systems matter: a single equation of order \(m\) becomes a system of \(m\) first-order equations, and then every
method above applies.
Note that the highest derivative never appears on the left of the
system — it is replaced by the equation, which is what makes
the reduction work.
Existence and Uniqueness
Definition 5.5
The function \(f(t,y_1,\ldots,y_m)\), defined on the set
\[
\begin{aligned}
D=\{(t,u_1,\ldots,u_m)\ |\ a\le t \le b, -\infty<u_i<\infty,
i=1,2,\ldots,m\}
\end{aligned}
\] is said to satisfy a Lipschitz condition on \(D\) in the variables \(u_1,u_2,\ldots,u_m\) if a constant \(L>0\) exists with \[
\begin{aligned}
|f(t,u_1,\ldots,u_m)-f(t,z_1,\ldots,z_m)| \le L\sum_{j=1}^m |u_j-z_j|,
\end{aligned}
\] for all \((t,u_1,\ldots,u_m)\) and \((t,z_1,\ldots,z_m)\) in \(D\).
Existence and Uniqueness
Theorem 5.8
Suppose \[
\begin{aligned}
D=\{(t,u_1,\ldots,u_m)\ |\ a\le t \le b, -\infty<u_i<\infty,
i=1,2,\ldots,m\}
\end{aligned}
\] and let \(f_i(t,u_1,\ldots,u_m)\), for each \(i=1,2,\ldots,m\), be continuous on \(D\) and satisfy a Lipschitz condition
there. The system of first-order differential equations \[
\begin{aligned}
\frac{du_k}{dt}(t) = f_k(t,u_1,\ldots,u_m),\quad u_k(a)=\alpha_k,\quad
k=1,\ldots,m
\end{aligned}
\] has a unique solution \(u_1(t),\ldots,u_m(t)\) for \(a\le t \le b\).
Numerical Methods
Numerical methods for systems of first-order differential equations
are vector-valued generalizations of methods for single equations.
Fourth order Runge-Kutta for systems
\[
\begin{aligned}
\mathbf{w}_0 &= \mathbf{\alpha} \\
\mathbf{k}_1 &= h \mathbf{f}(t_i,\mathbf{w}_i) \\
\mathbf{k}_2 &= h \mathbf{f}(t_i+\frac{h}{2},\mathbf{w}_i+\frac12
\mathbf{k}_1) \\
\mathbf{k}_3 &= h \mathbf{f}(t_i+\frac{h}{2},\mathbf{w}_i+\frac12
\mathbf{k}_2) \\
\mathbf{k}_4 &= h \mathbf{f}(t_{i+1},\mathbf{w}_i+\mathbf{k}_3) \\
\mathbf{w}_{i+1} &= \mathbf{w}_i + \frac16 (\mathbf{k}_1 +
2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4)
\end{aligned}
\] where \(\mathbf{w}_i=(w_{i,1},\ldots,w_{i,m})\) is
the vector of unknowns.
5.6 Multistep Methods
Multistep Methods
Definition 5.6
An \(m\)-step multistep
method for solving the initial-value problem \[
\begin{aligned}
y'=f(t,y),\quad a\le t \le b,\quad y(a)=\alpha,
\end{aligned}
\] has a difference equation for approximate \(w_{i+1}\) at \(t_{i+1}\): \[
\begin{aligned}
w_{i+1} &= a_{m-1}w_i + a_{m-2}w_{i-1} + \cdots + a_0w_{i+1-m} \\
&+ h[b_m f(t_{i+1},w_{i+1}) + b_{m-1} f(t_i,w_i) + \cdots \\
&+ b_0 f(t_{i+1-m},w_{i+1-m})],
\end{aligned}
\] where \(h=(b-a)/N\), and
starting values are specified: \[
\begin{aligned}
w_0=\alpha,\quad w_1=\alpha_1,\quad \ldots,\quad w_{m-1}=\alpha_{m-1}
\end{aligned}
\]Explicit method if \(b_m=0\), implicit method if \(b_m\ne 0\).
Runge-Kutta methods get their accuracy by evaluating \(f\) several times per step and then
throwing all of those values away. A multistep method
reuses the values it has already computed at previous nodes instead, so
a higher-order method can cost a single new evaluation per step.
Figure 5.9: Multistep
methods reuse past values of \(f\).
Above, the solution and the computed iterates; below, the slope
function, with the interpolating polynomial that AB2 and AB1 integrate
over the step. The iterates here are AB2’s own, so they drift slightly
from the solution.
It is second order but needs only one new evaluation of \(f\) per step, against two for a
second-order Runge-Kutta method. The price is that it needs \(w_1\) from somewhere before it can
start.
Integrate the initial-value problem \[
\begin{aligned}
y'=f(t,y),\quad a\le t \le b,\quad y(a)=\alpha
\end{aligned}
\] over \([t_i,t_{i+1}]\): \[
\begin{aligned}
y(t_{i+1}) = y(t_i) + \int_{t_i}^{t_{i+1}} f(t,y(t))\,dt
\end{aligned}
\] Replace \(f\) by polynomial
\(P(t)\) interpolating \((t_0,w_0),\ldots,(t_i,w_i)\), and
approximate \(y(t_i)\approx w_i\):
\[
\begin{aligned}
y(t_{i+1}) \approx w_i + \int_{t_i}^{t_{i+1}} P(t)\,dt
\end{aligned}
\]
The Adams methods come from integrating the differential equation
rather than differentiating it. Integrating \(y' = f(t,y)\) over one step, \[
y(t_{i+1}) = y(t_i) + \int_{t_i}^{t_{i+1}} f(t, y(t))\,dt,
\] which is exact. Now approximate the integrand by a polynomial
\(P(t)\) interpolating \(f\) at nodes we have already visited, and
replace \(y(t_i)\) by \(w_i\): \[
y(t_{i+1}) \approx w_i + \int_{t_i}^{t_{i+1}} P(t)\,dt .
\]
Every Adams method is this formula with a different choice of
interpolation nodes.
Figure 5.10: Left: AB2
interpolates the two most recent past values and integrates the
line beyond its data, so the method is explicit. Right: AM1 includes
\(f_{i+1}\), so the unknown appears on
both sides and the method is implicit.
Example 5.11
AB1 is Euler’s method.
Let \(P(t)\) interpolate the single
value \(f(t_i,w_i) \equiv f_i\), so
\(P\) is the constant \(f_i\). Then \[
\int_{t_i}^{t_{i+1}} P(t)\,dt = h f_i
\qquad\Longrightarrow\qquad
w_{i+1} = w_i + h\,f(t_i,w_i),
\] which is Euler’s method. So Euler is the one-step
Adams-Bashforth method.
Example 5.12
AB2, derived.
Let \(P(t)\) interpolate \(f_{i-1}\) and \(f_i\) — the two most recent values. In
Lagrange form, \[
P(t) = f_i\,\frac{t - t_{i-1}}{t_i - t_{i-1}}
+ f_{i-1}\,\frac{t - t_i}{t_{i-1} - t_i}
= f_i\,\frac{t - t_{i-1}}{h} - f_{i-1}\,\frac{t - t_i}{h}.
\]
Integrating from \(t_i\) to \(t_{i+1}\), \[
w_{i+1} = w_i + \left[f_i\frac{(t-t_{i-1})^2}{2h}
- f_{i-1}\frac{(t-t_i)^2}{2h}\right]_{t_i}^{t_{i+1}} .
\] At the upper limit \(t - t_{i-1} =
2h\) and \(t - t_i = h\); at the
lower limit the first bracket contributes \(h^2\) and the second vanishes. So \[
w_{i+1} = w_i + f_i\frac{(2h)^2}{2h} - f_{i-1}\frac{h^2}{2h}
- f_i\frac{h^2}{2h} + 0
= w_i + \frac{3h}{2}f(t_i,w_i) - \frac{h}{2}f(t_{i-1},w_{i-1}),
\] which is AB2. Note that the interpolation nodes \(t_{i-1}, t_i\) both lie outside
the interval of integration \([t_i,
t_{i+1}]\) — the polynomial is extrapolated, which is the source
of these methods’ step-size restrictions.
Derivation of Multistep Methods
Adams-Bashforth explicit \(m\)-step:
Newton backward-difference polynomial through \((t_i,f(t_i,y(t_i))),\ldots,(t_{i+1-m},f(t_{i+1-m},y(t_{i+1-m})))\).
\[
\begin{aligned}
\int_{t_i}^{t_{i+1}} f(t,y(t))\,dt &\approx
\int_{t_i}^{t_{i+1}} \sum_{k=0}^{m-1}
(-1)^k \binom{-s}{k} \nabla^k f(t_i,y(t_i))\,dt \\
&=\sum_{k=0}^{m-1} \nabla^k f(t_i,y(t_i)) h (-1)^k
\int_0^1 \binom{-s}{k}\,ds
\end{aligned}
\]
Now \(w_{i+1}\) appears on
both sides: this is an implicit method, and
each step requires solving an equation for \(w_{i+1}\) — by Newton’s method from Chapter
2, in general. Section 5.11 explains why that extra work is sometimes
unavoidable.
A one-step difference-equation with local truncation error \(\tau_i(h)\) is said to be
consistent if \[
\begin{aligned}
\lim_{h\rightarrow 0} \max_{1\le i \le N} |\tau_i(h)| = 0
\end{aligned}
\]
Definition 5.9
A one-step difference equation is said to be convergent if
\[
\begin{aligned}
\lim_{h\rightarrow 0} \max_{1\le i \le N} | w_i - y(t_i) | = 0,
\end{aligned}
\] where \(y_i=y(t_i)\) is the
exact solution and \(w_i\) the
approximation.
Two properties together give convergence. A method is
consistent if its local truncation error vanishes as
\(h \to 0\) — it solves the right
equation in the limit. It is stable if its own errors
do not amplify. Consistency alone is not enough, as we are about to
see.
Convergence of One-Step Methods
Theorem 5.9
Suppose the initial-value problem \(y'=f(t,y)\), \(a\le t \le b\), \(\quad
y(a)=\alpha\) is approximated by a one-step difference method in
the form \(w_0=\alpha\), \(w_{i+1}=w_i + h \phi(t_i,w_i,h)\). Suppose
also that \(h_0>0\) exists and \(\phi(t,w,h)\) is continuous with a
Lipschitz condition in \(w\) with
constant \(L\) on \(D\), then: \[
\begin{aligned}
D=\{(t,w,h)\ |\ a\le t \le b, -\infty<w<\infty, 0\le h\le h_0 \}.
\end{aligned}
\]
The method is stable;
The method is convergent if and only if it is consistent: \[
\begin{aligned}
\phi(t,y,0) = f(t,y)
\end{aligned}
\]
If \(\tau\) exists s.t. \(|\tau_i(h)|\le \tau(h)\) when \(0\le h \le h_0\), then \[
\begin{aligned}
|y(t_i) - w_i| \le \frac{\tau(h)}{L} e^{L(t_i-a)}.
\end{aligned}
\]
Proof of stability of the midpoint method.
The midpoint method is \(w_{i+1} = w_i +
h\phi(t_i, w_i, h)\) with \[
\phi(t,w,h) = f\!\left(t + \tfrac{h}{2},\; w +
\tfrac{h}{2}f(t,w)\right),
\] and we showed its local truncation error is \(\tau(h) = O(h^2)\), so it is consistent.
For stability we need \(\phi\) to be
Lipschitz in \(w\).
Suppose \(f\) is Lipschitz with
constant \(L\). Then \[
\begin{aligned}
|\phi(t,w,h) - \phi(t,\bar w,h)|
&= \left|f\!\left(t + \tfrac{h}{2}, w + \tfrac{h}{2}f(t,w)\right)
- f\!\left(t + \tfrac{h}{2}, \bar w + \tfrac{h}{2}f(t,\bar
w)\right)\right| \\
&\le L\left|w + \tfrac{h}{2}f(t,w) - \bar w - \tfrac{h}{2}f(t,\bar
w)\right| \\
&\le L\left[|w - \bar w| + \tfrac{h}{2}L|w - \bar w|\right] \\
&= \underbrace{\left(L + \tfrac{hL^2}{2}\right)}_{\textstyle
L'}|w - \bar w| ,
\end{aligned}
\] where the first inequality applies the Lipschitz condition in
the second argument and the second applies it again inside. So \(\phi\) is Lipschitz with constant \(L'\), and the method is stable.
For one-step methods consistency and stability together give
convergence, with \[
|y(t_i) - w_i| \le \frac{Kh^2}{L'}e^{L'(t_i - a)} .
\]■
Multistep methods can be
unstable
For multistep methods stability is a real constraint, not a
formality.
Example 5.14
A consistent method that blows up.
Consider \[
w_{n+1} = 3w_n - 2w_{n-1} - h\,f(t_{n-1}, w_{n-1}).
\] Writing \(y_n = y(t_n)\),
\(f_n = f(t_n, y(t_n))\) and \(y_n' = y'(t_n)\), the local
truncation error is \[
\tau_{n+1}(h) = \frac{y_{n+1} - 3y_n + 2y_{n-1}}{h} + f_{n-1} .
\]
So the method is consistent. Now watch it fail. Take
\(f \equiv 0\), for which the exact
solution is constant. The method becomes \[
w_{n+1} = 3w_n - 2w_{n-1}.
\] Looking for solutions \(w_n =
c\lambda^n\) with \(\lambda \ne
0\), \[
c\lambda^{n+1} = 3c\lambda^n - 2c\lambda^{n-1}
\qquad\Longrightarrow\qquad
\lambda^2 - 3\lambda + 2 = 0
\qquad\Longrightarrow\qquad
\lambda = 1,\ 2,
\] so the general solution is \[
w_n = c_1\cdot 1^n + c_2\cdot 2^n .
\] The \(2^n\) term means
any perturbation — a round-off error in \(w_1\), say — doubles at every step. The
method is consistent and useless.
Root Condition
Definition 5.10
Let \(\lambda_1,\ldots,\lambda_m\)
denote the roots of the characteristic equation \[
\begin{aligned}
P(\lambda) = \lambda^m - a_{m-1}\lambda^{m-1} - \cdots -
a_1\lambda - a_0 = 0
\end{aligned}
\] associated with the multistep method \[
\begin{aligned}
w_{i+1} =& a_{m-1}w_i + a_{m-2}w_{i-1} + \cdots + a_0 w_{i+1-m} \\
&+ hF(t_i,h,w_{i+1},w_i,\ldots,w_{i+1-m}).
\end{aligned}
\] If \(|\lambda_i|\le 1\) and
all roots with absolute value 1 are simple, the method is said to
satisfy the root condition.
The example shows what to check: the roots of the
characteristic polynomial must not exceed \(1\) in modulus. That is the root
condition.
Example 5.15
AB3 satisfies the root condition.
For \[
w_{i+1} = w_i + \frac{h}{12}\left(23f_i - 16f_{i-1} + 5f_{i-2}\right),
\] an \(m = 3\) step method, one
can show \(\tau(h) = O(h^3)\), so it is
consistent. Its characteristic polynomial is \[
P(\lambda) = \lambda^m - \lambda^{m-1} = \lambda^3 - \lambda^2
= \lambda^2(\lambda - 1),
\] with roots \[
\lambda = 0 \text{ (multiplicity 2)}, \qquad \lambda = 1 \text{
(multiplicity 1)}.
\] Only \(\lambda = 1\) has
modulus \(1\) and it is simple, while
the others are strictly inside the unit circle — so AB3 satisfies the
root condition and is strongly stable. Consistency plus
stability gives convergence.
Stability
Definition 5.11
Methods that satisfy the root condition and have \(\lambda=1\) as the only root of magnitude
one are called strongly stable.
Methods that satisfy the root condition and have more than one
distinct root with magnitude one are called weakly stable.
Methods that do not satisfy the root condition are
unstable.
Theorem 5.10
A multistep method \[
\begin{aligned}
w_{i+1}=&a_{m-1}w_i + a_{m-2}w_{i-1} + \cdots + a_0 w_{i+1-m} \\
&+h F(t_i,h,w_{i+1},w_i,\ldots,w_{i+1-m})
\end{aligned}
\] is stable if and and only if it satisfies the root condition.
If it is also consistent, then it is stable if and only if it is
convergent.
5.11 Stiff Differential
Equations
Stiff Equations
A stiff differential equation is numerically unstable
unless the step size is extremely small
Large derivatives give error terms that are dominating the
solution
Example: The initial-value problem \[
\begin{aligned}
y'=-30y, \quad 0\le t \le 1.5, \quad y(0)=\frac13
\end{aligned}
\] has exact solution \(y=\frac13
e^{-30t}\). But RK4 is unstable with step size \(h=0.1\).
Some problems restrict the step size for reasons that have nothing to
do with accuracy.
Example 5.16
Euler on \(y' = -y\) with too
large a step.
The exact solution of \(y' =
-y\), \(y(0) = 1\) is \(y(t) = e^{-t}\), which decays
monotonically. Take Euler with \(h =
2\) on \(0 \le t \le 10\): \[
\begin{aligned}
w_0 &= 1, \\
w_1 &= w_0 + h f(t_0,w_0) = w_0 + h(-w_0) = (1-h)w_0 = -1, \\
w_2 &= (1-h)w_1 = +1,
\end{aligned}
\] and in general \[
w_i = (1-h)^i w_0 .
\]
Figure 5.11: Euler’s method
on \(y'=-y\) with \(h=2\). The exact solution decays; the
numerical one oscillates between \(\pm
1\) and never settles.
With \(h = 2\) the factor \(1 - h = -1\), so the iterates alternate
between \(+1\) and \(-1\) forever. The solution they are
supposed to approximate has decayed to \(e^{-10} \approx 5\times 10^{-5}\). Nothing
is wrong with the accuracy of Euler’s method here; the step
size is simply outside the range where the method is stable.
The test equation
Euler’s Method for Test Equation
Consider the simple test equation \[
\begin{aligned}
y'=\lambda y,\quad y(0)=\alpha, \quad \text{where } \lambda<0
\end{aligned}
\] with solution \(y(t)=\alpha
e^{\lambda t}\).
Stability requires \(|1+h\lambda|<1\), or \(h<2/|\lambda|\).
To make this precise, examine the test problem\[
y' = \lambda y, \qquad y(0) = \alpha, \qquad
\operatorname{Re}\lambda < 0,
\] whose exact solution \(\alpha
e^{\lambda t}\) decays. Euler’s method gives \[
w_i = \underbrace{(1 + h\lambda)}_{\textstyle
Q(h\lambda)}{}^{\,i}\,\alpha,
\] so the numerical solution decays if and only if \[
|1 + h\lambda| < 1 .
\] For real \(\lambda < 0\)
this is \(h < 2/|\lambda|\) —
precisely the restriction the example ran into, since \(\lambda = -1\) there gives \(h < 2\).
Regions of absolute
stability
Region of Stability
Definition 5.12
The region \(R\) of absolute
stability for a one-step method is \(R=\{
h\lambda\in\mathcal{C}\ | \ |Q(h\lambda)|<1\}\), and for a
multistep method, it is \(R=
\{h\lambda\in\mathcal{C}\ |\ |\beta_k|<1,\text{ for all
zeros }\beta_k\text{ of }Q(z,h\lambda)\}\).
A numerical method is said to be A-stable if its region
\(R\) of absolute stability contains
the entire left half-plane.
The set of \(h\lambda\) for which
\(|Q(h\lambda)| < 1\) is the
method’s region of absolute stability. For Euler it is
the disk of radius \(1\) centred at
\(-1\).
Example 5.17
The region for RK4.
Substituting \(f = \lambda y\) into
the four stages of RK4, \[
\begin{aligned}
k_1 &= h\lambda w_i, \\
k_2 &= h\lambda\!\left(w_i + \tfrac{k_1}{2}\right)
= h\lambda\!\left(w_i + \tfrac{h\lambda}{2}w_i\right)
= \left[h\lambda + \tfrac{(h\lambda)^2}{2}\right]w_i,
\end{aligned}
\] and continuing through \(k_3\) and \(k_4\), \[
w_{i+1} = \underbrace{\left[1 + h\lambda + \frac{(h\lambda)^2}{2}
+ \frac{(h\lambda)^3}{6} + \frac{(h\lambda)^4}{24}\right]}_{\textstyle
Q(h\lambda)}w_i .
\]
The bracket is the degree-four truncation of \(e^{h\lambda}\), which is no accident: a
method of order \(p\) must match the
exponential to order \(p\).
Figure 5.12: Regions of
absolute stability. Euler’s is the disk \(|1+h\lambda|<1\); RK4’s is the larger
lobe where the degree-four truncation of the exponential has modulus
less than one.
RK4’s region extends to about \(h\lambda =
-2.78\) on the negative real axis against Euler’s \(-2\), so RK4 buys a little more stability
along with its extra order — but only a little. That is the point of a
stiff problem: when \(\lambda\) is very large and negative, no
explicit method has a useful step size, and an implicit method such as
the trapezoidal method of §5.6 is required.
Multistep Methods
Apply a multistep method to the test equation: \[
\begin{aligned}
w_{j+1} =& a_{m-1}w_j + \cdots + a_0 w_{j+1-m} \\
&+h\lambda(b_mw_{j+1} + b_{m-1} w_j + \cdots
+ b_0 w_{j+1-m})
\end{aligned}
\] or \[
\begin{aligned}
(1-h\lambda b_m)w_{j+1} - (a_{m-1}+h\lambda b_{m-1}) w_j
- \cdots - (a_0 + h\lambda b_0) w_{j+1-m} = 0
\end{aligned}
\] Let \(\beta_1,\ldots,\beta_m\) be the zeros of
the characteristic polynomial\[
\begin{aligned}
Q(z,h\lambda) = (1-h\lambda b_m)z^m - (a_{m-1}+h\lambda b_{m-1}) z^{m-1}
-\cdots - (a_0+h\lambda b_0)
\end{aligned}
\] Then \(c_1,\ldots,c_m\) exist
with \[
\begin{aligned}
w_j = \sum_{k=1}^m c_k (\beta_k)^j
\end{aligned}
\] and \(|\beta_k|<1\) is
required for stability.
