Lecture 11 · Ordinary Differential Equations

COMP 521 · Lecture 11

Ordinary differential equations

Models, solutions, and numerical methods

Anna Rosen · San Diego State University

What an ODE describes

What does an ODE tell us?

An ODE relates an unknown function to its derivatives with respect to one independent variable.

\frac{dy}{dx}=f(x,y)

Independent variable

x

Time, depth, radius, …

Unknown state function

y(x)

State at coordinate $x$

Derivative

\frac{dy}{dx}

Change per unit coordinate

Think: Does a slope rule determine a unique function? What else is needed?

What an ODE describes

What information belongs in the state?

The state is the set of quantities we track to describe a system.

It can be a scalar $y$ or a vector $\mathbf y$ of coupled variables.

\frac{d\mathbf y}{dx}=\mathbf f(x,\mathbf y)

Examples:

Population · scalar

y(t)=N(t)

How many individuals?

Spring · vector

A mass on a spring; displacement is measured from equilibrium.

\mathbf y(t)=\begin{pmatrix}x(t)\\v(t)\end{pmatrix}

Position and velocity

Star or planet · vector

\mathbf y(r)=\begin{pmatrix}m(r)\\P(r)\end{pmatrix}

Enclosed mass and pressure

Predict: Two spring masses have the same position but opposite velocities. How will their positions change next?

From a system to a mathematical model

Modeling: turn assumptions into equations

Choose the quantities. What do we want to describe?
State the assumptions. What detail can we neglect?
Write relationships. What laws or rules connect the quantities?
Set parameters and conditions. What specifies the problem?

From assumptions to a rate law

An ODE model: radioactive decay

Example: one radioactive isotope

State: $N(t)$ , expected undecayed nuclei.
Assumptions: independent decays; no replenishment.

\text{Rate law:}\qquad \frac{dN}{dt}=-\lambda N

Parameter: constant $\lambda>0$ , units $\mathrm{s^{-1}}$ .

Initial condition: $N(0)=N_0$ .

Initial value problem (IVP): an ODE and its starting state.

ODE models describe change with one independent variable, such as time or radius.

Applications: reactions and growth in time; pressure balance along radius.

What an ODE describes

“Ordinary” and “first order” mean different things

Ordinary: derivatives with respect to one independent variable.
Order: the highest derivative present.

First order

\frac{dN}{dt}=rN

Population growth

Second order

m\frac{d^2x}{dt^2}=-kx

An undamped spring: mass $m$ , stiffness $k$ (force per displacement).

Both are ODEs. Neither “ordinary” nor “order” requires time as the coordinate.

Define an ODE and its state

Partial differential equations (PDEs)

Heat equation

\frac{\partial T}{\partial t}=\alpha\frac{\partial^2T}{\partial x^2}

$T(x,t)$ : temperature varies with position and time.
$\alpha$ : thermal diffusivity, in $\mathrm{cm^2\,s^{-1}}$ .

One-temperature ODE: $T(t)$ keeps one value for the whole object.

Temperature-excess profiles at scaled times 0, 0.5, and 1 on an infinite rod. The peak becomes lower and the hot region becomes wider as heat diffuses. — Heat diffusion lowers and broadens the peak. Scaled variables.

PDEs are central to computational modeling: fluid dynamics, electrostatics, and general relativity.

From a rate law to an exact solution

Model a growing bacterial culture

State: $N(t)$ , population size.

Rate law:

\frac{dN}{dt}=rN

Parameter: constant net per-cell rate $r>0$ ( $\mathrm{s^{-1}}$ ).
Initial condition: $N(0)=N_0>0$ .

Assumptions: no resource limit or migration.

Exact normalized population growth for a constant per-capita rate; rt is dimensionless. — Early culture growth with a constant net rate per cell.

Computed exponential population growth with dimensionless rt and normalized population axes. — Precomputed Python result

Exact growth and decay models

Derive exponential growth and doubling time

$r$ : constant net growth rate per individual, units $\mathrm{time^{-1}}$ .

For $N(0)=N_0>0$ and $r>0$ :

\frac{dN}{N}=r\,dt

\begin{aligned} \int_{N_0}^{N(t)}\frac{dN}{N}&=\int_0^t r\,d\tau\\ \ln\frac{N(t)}{N_0}&=rt \end{aligned}

N(t)=N_0e^{rt}

Doubling time: set $N(t_2)=2N_0$ .

2=e^{rt_2}\quad\Rightarrow\quad t_2=\frac{\ln2}{r}

Exponential growth doubles at each equal interval of length ln(2)/r. — Equal doubling intervals add equal factors, not equal numbers.

Exact growth and decay models

Derive exponential decay and half-life

Rate law: $N'=-\lambda N$ , with constant $\lambda>0$ and $N(0)=N_0>0$ .

Separate:

\frac{dN}{N}=-\lambda\,dt

\begin{aligned} \int_{N_0}^{N(t)}\frac{dN}{N}&=-\lambda\int_0^t d\tau\\ \ln\frac{N(t)}{N_0}&=-\lambda t \end{aligned}

N(t)=N_0e^{-\lambda t}

Half-life: set $N(t_{1/2})=N_0/2$ .

\ln(1/2)=-\lambda t_{1/2}\quad\Rightarrow\quad t_{1/2}=\frac{\ln2}{\lambda}

Remaining population on a logarithmic y axis. Equal half-life intervals halve the population and form a straight line. — Log y: equal time intervals give equal fractional drops.

Growth and decay on log y: straight lines with opposite slopes, starting from the same normalized population.

$a=r>0$ : growth. $a=-\lambda<0$ : decay.

Exact growth and decay models

Derive cooling toward the bath temperature

Bath: fixed $T_{\rm cold}$ .
Initial state: $T_{\rm hot}(0)=T_0>T_{\rm cold}$ .

1. Shift the temperature

\begin{aligned} \theta(t)&=T_{\rm hot}(t)-T_{\rm cold}\\ \theta_0&=T_0-T_{\rm cold} \end{aligned}

2. Separate and integrate

\begin{aligned} \frac{d\theta}{\theta}&=-k\,dt\\ \ln\frac{\theta(t)}{\theta_0}&=-kt \end{aligned}

3. Recover the temperature

\begin{aligned} \theta(t)&=\theta_0e^{-kt}\\ T_{\rm hot}(t)&=T_{\rm cold}+\theta_0e^{-kt} \end{aligned}

Limit: $\theta\to0$ ; $T_{\rm hot}\to T_{\rm cold}$ .

Newton’s law of cooling

\frac{dT_{\rm hot}}{dt}=-k(T_{\rm hot}-T_{\rm cold})

Uniform temperature within the object.
Fixed bath temperature $T_{\rm cold}$ .
Constant rate: $k>0$ , units $\mathrm{s^{-1}}$ .

Logarithmic normalized temperature excess. Exponential cooling is a straight line and remains positive while approaching equilibrium. — Log y reveals the cooling tail; the temperature excess stays positive.

Spatial ODEs and computational motivation

Derive pressure as a function of depth

\frac{dP}{dz}=\rho g

$z$ increases downward from the surface.
Assume constant $\rho$ and $g$ ; $P(0)=P_{\rm surface}$ .

\begin{aligned} \int_{P_{\rm surface}}^{P(z)} dP&=\int_0^z\rho g\,ds\\ P(z)-P_{\rm surface}&=\rho gz\\ P(z)&=P_{\rm surface}+\rho gz \end{aligned}

Explain: Why can $dP/dz\ne0$ while $\partial P/\partial t=0$ ?

Parameters: $\rho$ : density; $g$ : gravitational acceleration.

Pressure increases linearly with downward depth when density and gravity are constant. — Pressure rises linearly with depth for constant density and gravity.

Spatial ODEs and computational motivation

Inside a planet, radius is the coordinate

\mathbf y(r)=\begin{pmatrix}m(r)\\P(r)\end{pmatrix}

\begin{aligned}\frac{dm}{dr}&=4\pi r^2\rho\\\frac{dP}{dr}&=-\frac{Gm\rho}{r^2}\end{aligned}

$m(r)$ : enclosed mass.
$P(r)$ : pressure.
$G$ : gravitational constant.

Outward: mass increases; pressure decreases.

A highlighted local patch of a spherical shell, with inner pressure force outward, outer pressure force inward, and gravity inward.

Hydrostatic equilibrium: pressure forces balance gravity.

Newton’s second law: $F_{\rm pressure}+F_{\rm gravity}=0$ .

Static model: neither expanding nor contracting.

Spatial ODEs and computational motivation

The equations need a material model and boundaries

\rho=\rho(P)

m(0)=0,\qquad P(R)=0

Equation of state: pressure determines density, $\rho(P)$ .

Start with a trial central pressure; find the surface $R$ .

Diagnose: Can smaller radial steps fix an incorrect $\rho(P)$ ? Explain.

Density closure joins radial mass and pressure equations into a specified planet model.

Spatial ODEs and computational motivation

We know the slope. Why still compute?

\frac{dy}{dx}=f(x,y)

Knowing the derivative does not guarantee a useful solution formula.
Coupled equations, nonlinear material laws, and events lead us to compute.

Computable derivatives lead to an approximate solution with explicit numerical checks.

Spatial ODEs and computational motivation

Specify an initial value problem (IVP)

IVP: an ODE plus a specified starting state.

\begin{aligned} \frac{d\mathbf y}{dx}&=\mathbf f(x,\mathbf y) &&\text{ODE}\\ \mathbf y(x_0)&=\mathbf y_0 &&\text{initial condition} \end{aligned}

$x_0$ : starting coordinate.
$\mathbf y_0$ : full starting state.
Solve: find $\mathbf y(x)$ satisfying both over the chosen interval.

Predict

For $y'=y$ , if $y(0)$ doubles, does $y(1)$ double too? Why?

Explain

Yes: $2y(x)$ obeys the same ODE with twice the starting value.

Two initial values select different exact solutions in the same slope field dy/dx=y. — Two initial values select different solutions of the same ODE.

From slopes to numerical steps

Approximate a derivative on a grid

Grid: $x_n=x_0+nh$ , with $h>0$ . Assume smooth $y(x)$ and exact samples.

Forward: derivative error $O(h)$ .

\frac{y(x_n+h)-y(x_n)}{h}

Backward: derivative error $O(h)$ .

\frac{y(x_n)-y(x_n-h)}{h}

Centered: derivative error $O(h^2)$ .

\frac{y(x_n+h)-y(x_n-h)}{2h}

Schematic of three equally spaced samples on a smooth curve. Forward uses the middle and right samples, backward uses the left and middle, and centered uses the two outer samples. — Nearby values give secant slopes that approximate $y'(x_n)$ .

Euler, stability, and error

Euler turns a slope into a step

Euler: first order. Global state error $O(h)$ over a fixed interval.

Discretize the derivative:

\frac{y_{n+1}-y_n}{h}=f(x_n,y_n)

$y_n\approx y(x_n)$ : a computed state.

Solve for the next state:

y_{n+1}=y_n+h f(x_n,y_n)

Use the starting slope for the whole step.

Predict

For growth $y'=y$ , $y(0)=1$ , the exact solution is $y(t)=e^t$ . Will Euler lie above or below this curve? Why?

Decay model: $u'=-u$ , $u(0)=1$ .

Exact solution: $u_{\mathrm{exact}}(t)=e^{-t}$ .

Euler follows the initial tangent for u′=−u; the convex exact decay lies above it.

Exact increment: integrate the ODE

\begin{aligned} \Delta y&=y(x_n+h)-y(x_n)\\ &=\int_{x_n}^{x_n+h}f(s,y(s))\,ds \end{aligned}

The intermediate state $y(s)$ is unknown.

Euler uses a left rectangle:

\Delta y\approx h f(x_n,y_n)

Starting slope × step approximates the increment.

Euler, stability, and error

Predict one Euler step

Decay model and exact solution:

\begin{aligned} u'&=-u,\qquad u(0)=u_0=1\\ u_{\mathrm{exact}}(t)&=e^{-t} \end{aligned}

Step: $t_0=0$ , $h=0.2$ , so $t_1=0.2$ .

Predict

Compute the Euler value $u_1$ . Compare it with $u_{\mathrm{exact}}(0.2)$ from the displayed solution. Explain using curvature.

\begin{aligned} f(0,u_0)&=-u_0=-1\\ u_1&=1+0.2(-1)=0.8 \end{aligned}

Exact: $u_{\mathrm{exact}}(0.2)=e^{-0.2}\approx0.818731$ .

One Euler step from u(0)=1 with h=0.2 gives 0.8; the exact endpoint is about 0.818731. — The starting tangent undershoots the convex exact decay.

Euler, stability, and error

Shrinking magnitude does not guarantee positive states

Same decay model: $u'=-u$ , $u(0)=1$ .

Exact: $u_{\mathrm{exact}}(t)=e^{-t}>0$ .

Euler: substitute $f(t_n,u_n)=-u_n$ .

u_{n+1}=u_n-hu_n=(1-h)u_n

$h=0.25$ : multiplier 0.75 → positive decay.

$h=1.5$ : -0.5 → alternating sign. $h=3$ : -2 → growing magnitude.

Euler with h=1.5 alternates sign while its magnitude shrinks; exact decay stays positive. — Linear axes retain the sign changes.

Log magnitude of exact decay and Euler at steps 0.25, 1.5 and 3 over the same time interval. Step 3 grows while the smaller steps shrink at different rates; this panel hides sign. — Log magnitude separates decay rates and instability. Signs are shown in the preceding plot.

Euler, stability, and error

One old state supplies every coupled slope

\mathbf y_{n+1}=\mathbf y_n+h\,\mathbf f(x_n,\mathbf y_n)

Compute every slope from the same old state.

Slope → increment → new state.

Euler computes all coupled slopes from the same old state.

Euler, stability, and error

What does Euler leave out?

Exact step: Taylor expansion

\begin{aligned} y(x+h)&=y(x)+h y'(x)\\ &\quad+\tfrac12 h^2 y''(x)+O(h^3) \end{aligned}

The ODE supplies $y'=f(x,y)$ .

Euler keeps the constant and linear terms:

y_{\rm Euler}=y(x)+h f(x,y(x))

First omitted term:

\tfrac12 h^2y''(x)

One-step state defect: $O(h^2)$ .
Divide by $h$ : derivative defect $O(h)$ .

Over a fixed interval: global error $O(h)$ .

Decay model: $u'=-u$ , $u(0)=1$ .

Exact solution: $u_{\mathrm{exact}}(t)=e^{-t}$ .

Euler, stability, and error

Local defects become global error

State error: $e_n=\|\mathbf y_n-\mathbf y(x_n)\|$ .

$L$ : bound on the slope’s sensitivity to the state.
$Ch^2$ : bound on the new defect from one Euler step.

e_{n+1}\leq(1+hL)e_n+Ch^2

\begin{aligned} e_0&=0,\quad T=Nh\ \mathrm{fixed}\\ e_N&=O(h) \end{aligned}

Propagated error + new defects: first-order global error.

Decay model: $u'=-u$ , $u(0)=1$ .

Exact solution: $u_{\mathrm{exact}}(t)=e^{-t}$ .

Five Euler steps and exact decay on focused linear axes. A vertical arrow marks the endpoint state error at T=1.

Euler, stability, and error

Halve the step; measure the error

Decay benchmark: $u'=-u$ , $u(0)=1$ .

\begin{aligned} u_{\mathrm{exact}}(t)&=e^{-t}\\ e(h)&=|u_h(T)-e^{-T}| \end{aligned}

\frac{e(h)}{e(h/2)}\approx2^p

$p$ : the method’s global order.

Hold fixed: initial state, endpoint $T$ , and error norm.

Euler: halving $h$ should approach half the error.

Coarse and fine grids end at exactly the same coordinate; compare states at that endpoint. — Both grids end at the same coordinate.

Oscillators as numerical laboratories

A spring needs position and velocity

m\frac{d^2x}{dt^2}=-kx,\qquad\omega=\sqrt{k/m}

\frac{dx}{dt}=v,\qquad\frac{dv}{dt}=-\omega^2x

$m$ : mass; $k$ : spring stiffness, force per displacement.
$x$ : displacement from equilibrium; $v$ : velocity.
$\omega$ : angular frequency, units $\mathrm{s^{-1}}$ .

Position and velocity determine the next change.

A mass m on a spring of stiffness k is displaced right from equilibrium x=0; the restoring force F=-kx points left. — **Hooke’s law:** the restoring force opposes displacement.

The harmonic oscillator as a benchmark

Derive the oscillator benchmark

Displacement amplitude: $A$ .

Dimensionless variables: $q=\frac{x}{A}$ , $p=\frac{v}{A\omega}$ , $t=\omega t_{\rm phys}$ .

q'=p,\qquad p'=-q

Initial state: $q(0)=1$ , $p(0)=0$ .

\begin{aligned} q''+q&=0\\ q(t)&=C\cos t+D\sin t \end{aligned}

\begin{aligned}C&=q(0)=1\\D&=q'(0)=p(0)=0\end{aligned}

\begin{aligned}q(t)&=\cos t\\p(t)&=q'(t)=-\sin t\end{aligned}

Exact cosine position and negative-sine velocity for the scaled oscillator with initial state (1,0). — A quarter-cycle phase offset.

The harmonic oscillator as a benchmark

Phase diagram: same energy, different states

Phase diagram: each point $(q,p)$ represents a state.
Energy conservation imposes an algebraic constraint:

\begin{aligned} E=\tfrac12(q^2+p^2)&=\tfrac12\\ q^2+p^2&=1 \end{aligned}

$(1,0)$ and $(0,-1)$ lie on the same energy circle.

Energy fixes the circle, not the state at a given time.

Phase diagram with scaled position q on the horizontal axis and scaled velocity p on the vertical axis. Distinct states lie on the same constant-energy circle. — Horizontal: scaled position. Vertical: scaled velocity. Equal axis scales.

The harmonic oscillator as a benchmark

Why does Euler spiral outward?

\begin{aligned}q_{n+1}&=q_n+hp_n\\p_{n+1}&=p_n-hq_n\end{aligned}

E_{n+1}=(1+h^2)E_n

Both updates use the original state.
Phase experiment: $(q_0,p_0)=(1,0)$ , $h=0.1$ , $T=12$ .

Every $h\ne0$ expands this numerical orbit.

Equal-scale q,p phase plane. Euler h=0.1 over scaled duration 12 spirals outward. — Euler spirals outward; the exact oscillator stays on the circle.

Midpoint and classical RK4

Midpoint: predict, sample, update

Explicit midpoint (RK2): second order. Global state error $O(h^2)$ .

\begin{aligned}k_1&=f(x_n,y_n)\\k_2&=f(x_n+h/2,y_n+hk_1/2)\end{aligned}

y_{n+1}=y_n+hk_2

Predict a temporary midpoint state.
Update from the original state, using the midpoint slope.

Predict

For $y'=y$ , set $k_2=f(x_n+h/2,y_n)$ . Does the update differ from Euler?

For $y'=y$ , missing the midpoint state gives Euler, first order.

Same oscillator: $q'=p$ , $p'=-q$ .

$(q_0,p_0)=(1,0)$ ; $h=0.1$ , $T=12$ .
RK2: much less outward drift than Euler.

E_{n+1}=\left(1+\frac{h^4}{4}\right)E_n

Small drift is still not exact energy conservation.

Explicit midpoint: starting slope, predicted midpoint, new slope, and old-state base.

Local Taylor check: set $y_n=y(x_n)=y$ .

\begin{aligned} y''&=f_x+f_y f\\ k_2&=f+\tfrac h2 y''+O(h^2)\\ y_{n+1}&=y+h f+\tfrac{h^2}{2}y''+O(h^3) \end{aligned}

$f_x=\partial f/\partial x$ , $f_y=\partial f/\partial y$ .
Matches the exact step through $h^2$ .
One-step $O(h^3)$ → global $O(h^2)$ .

Phase diagram cropped to plus or minus 1.08 on both axes. RK2 is bold teal; exact is dark dotted; Euler is faint gray and continues beyond this crop. A square marks the RK2 endpoint. — Near-orbit view; Euler continues outside this crop. Next: endpoint zoom.

Equal-scale endpoint zoom at T=12 separates the teal RK2 square from the dark exact point and orbit. The experiment is unchanged. — Endpoint zoom resolves the RK2 offset from exact at the same time.

Midpoint and classical RK4

Trace one midpoint step

y^\prime=y,\quad y_0=1,\quad h=0.1

Predict

Compute $y_{\rm mid}$ and $y_1$ . Is the trial midpoint exact at $x_0+h/2$ ?

k_1=1,\quad y_{\mathrm{mid}}=1.05,\quad k_2=1.05

Midpoint: $y_1=1+0.1(1.05)=1.105$ .
Exact midpoint: $e^{0.05}\approx1.051271096$ .
Exact endpoint: $e^{0.1}\approx1.105170918$ .

Log absolute state differences at the trial midpoint x=0.05 and final endpoint x=0.1. Trial error is about 0.00127; endpoint error is about 0.000171. — Log y resolves both approximation errors at their own coordinates.

Midpoint and classical RK4

RK4: build four trial slopes

Classical RK4: fourth order. Global state error $O(h^4)$ .

Trace: $y'=y$ , $y_0=1$ , $h=0.1$ .

k_1=f(x_n,y_n)

k_2=f(x_n+h/2,y_n+hk_1/2)

Predict

If we stop here and use $y_1=y_0+h k_2$ , what method and order do we get?

$k_1=1$ , $k_2=1.05$ .

Stopping here gives midpoint RK2, second order.

Known: $k_2=1.05$ .

k_3=f(x_n+h/2,y_n+hk_2/2)

Predict

For this same step, is $k_3=k_2$ ? Explain before calculating.

k_4=f(x_n+h,y_n+hk_3)

$k_3=1.0525$ , $k_4=1.10525$ .

Same midpoint coordinate; different predicted states.

RK4's starting slope supplies the first midpoint state; the displayed trace uses h=0.1. — The first trial midpoint.

RK4 stages 3 and 4 construct new trial states from earlier slopes, retaining the original state as their base. — New slopes use new trial states.

Midpoint and classical RK4

RK4: combine slopes, then update once

\bar k=\tfrac16(k_1+2k_2+2k_3+k_4)

y_{n+1}=y_n+h\bar k

\begin{aligned}y_1&=1+\tfrac{0.1}{6}(6.31025)\\&=1.105170833\ldots\end{aligned}

Exact midpoint: $e^{0.05}\approx1.051271096$ .
Exact endpoint: $e^{0.1}\approx1.105170918$ .
Update the original state once.

For $y'=y$ , the RK4 update becomes

\begin{aligned} R_{\rm RK4}(h)&=1+h+\tfrac{h^2}{2}\\ &\quad+\tfrac{h^3}{6}+\tfrac{h^4}{24}\\ y_{n+1}&=y_n R_{\rm RK4}(h) \end{aligned}

$R_{\rm RK4}$ : the step multiplier for this model.

Same oscillator: $q'=p$ , $p'=-q$ .

$(q_0,p_0)=(1,0)$ ; $h=0.1$ , $T=12$ .
RK4: nearly overlaps the exact circle.

E_{n+1}=\left(1-\frac{h^6}{72}+\frac{h^8}{576}\right)E_n

Energy is not conserved exactly.
Inspect the endpoint for state and phase error.

Completed scalar RK4 trace with the four stage states and slopes.

Exact step: $y(x+h)=y(x)e^h$ .

\begin{aligned} e^h&=1+h+\tfrac{h^2}{2}\\ &\quad+\tfrac{h^3}{6}+\tfrac{h^4}{24}+O(h^5) \end{aligned}

Matches the exact step through $h^4$ .
One-step state defect: $O(h^5)$ .
Global state error: $O(h^4)$ .

Order: global accuracy, not ODE order or stage count. This example is not a general RK4 order proof.

Phase diagram cropped to plus or minus 1.08 on both axes. Bold teal RK2 and purple RK4 use interleaved dashes and distinct endpoint markers. Euler is faint gray and continues beyond the crop; exact is dark dotted. — Teal RK2; purple RK4. Euler is faded and cropped. Next: endpoint zoom.

Magnified phase-plane endpoint at T=12 on equal axis scales. RK2 has a visible state offset; RK4 and exact nearly overlap at this magnification. — Endpoint zoom: RK2 is displaced. RK4 and exact nearly overlap at this scale.

Further phase-plane zoom centered on the exact endpoint at T=12. Equally scaled coordinate offsets in microunits separate RK4 from exact; an arrow shows the small state error. — Further zoom: offsets from the exact endpoint in millionths. RK4 has a small, nonzero state error.

Midpoint and classical RK4

Higher order versus more work per step

e=\|\mathbf y_h(2)-\mathbf y_{\mathrm{exact}}(2)\|_\infty

Same experiment: $(q_0,p_0)=(1,0)$ , $T=2$ ; halve $h$ .

Euler / midpoint / RK4: orders $1,2,4$ ; calls/step $1,2,4$ .
Compare work at equal error.

Common-endpoint oscillator errors at T=2 using h=0.1 halved four times; linear reference maps only. — State error at one common endpoint decreases with step refinement.

Two ways to improve a numerical answer

Richardson: refine, then cancel the leading error

Same problem, same endpoint: $u'=-u$ , $u_0=1$ , $h=0.2$ .

One full Euler step, or two half steps.

Predict

Are two half Euler steps a second-order method? Distinguish refinement from correction.

Y_{\mathrm{coarse}}=1-0.2=0.8

Y_{\mathrm{fine}}=(1-0.1)^2=0.81

Two half steps still give a first-order method.

Y_{\mathrm{ext}}=Y_{\mathrm{fine}}+\frac{Y_{\mathrm{fine}}-Y_{\mathrm{coarse}}}{2^p-1}

p=1:\quad Y_{\mathrm{ext}}=2(0.81)-0.8=0.82

This Euler correction is second order.

Assume a smooth solution and stable error propagation.
Order depends on the base method and its error expansion.

Euler step doubling over duration 0.2 gives coarse 0.8, fine 0.81 and corrected 0.82, compared with exact exp(-0.2). — Refinement and correction are different operations at the same endpoint.

Leapfrog, energy, and phase

Leapfrog: kick, drift, kick

Leapfrog belongs to the symplectic family: second order, global state error $O(h^2)$ .

p_{n+1/2}=p_n+\frac h2F(q_n)

$H=p^2/(2m)+V(q)$ ; force $F(q)=-V^\prime(q)$ .

q_{n+1}=q_n+\frac hm p_{n+1/2}

p_{n+1}=p_{n+1/2}+\frac h2F(q_{n+1})

State: $q$ : position; $p$ : momentum.

Energy: $H$ : total; $V(q)$ : potential; $m$ : mass.

Scaled oscillator: $m=1$ , $F(q)=-q$ .

The symmetric composition changes momentum, position, then momentum.

Leapfrog, energy, and phase

Preserve the structure of the equations

$\Phi_h$ : the one-step map from $(q_n,p_n)$ to $(q_{n+1},p_{n+1})$ .
$J=D\Phi_h$ : the state Jacobian of the map.

\Omega=\begin{pmatrix}0&1\\-1&0\end{pmatrix},\qquad J^{\mathsf T}\Omega J=\Omega

For one canonical pair: oriented phase-space area is preserved.

Symplectic structure, state accuracy, and physical energy are separate claims.

Fixed-step kick–drift–kick: symplectic, second order, and time reversible.

Leapfrog, energy, and phase

Same oscillator, different energy behavior

Resolution: 100 steps per period.
Duration: 20 periods.
Initial state: $(q_0,p_0)=(1,0)$ .
Same step; costs differ.

First plot (log): $|E-E_0|/E_0$ . Zero errors omitted.

RK2: growing error.
Leapfrog: bounded oscillations.
RK4: much smaller drift.

Next: signed RK4 error.

Logarithmic absolute relative energy error over 20 periods. RK2 grows, leapfrog oscillates below a bounded envelope, and RK4 has much smaller nonzero drift. Zero errors are omitted. — Energy-error magnitude on a logarithmic scale. Zero values are omitted; same step, unequal derivative cost.

Magnified nonzero RK4 energy error over 20 oscillator periods. — RK4 has a small, nonzero energy drift. Magnified vertical scale.

Leapfrog, energy, and phase

Good energy; wrong timing?

Phase: timing along the orbit.
Energy: its level set.
Log scale: phase-error magnitude. Zero values are omitted.
Same 100 steps per period and 20 periods.

Check the full state at the same time.

Logarithmic phase-error magnitude. RK2 and leapfrog accumulate larger timing errors than RK4 in the same oscillator experiment. Zero errors are omitted. — Energy alone misses accumulating timing error.

Return to Project 4 and synthesize

A planet model wraps integration in a root solve

P(R)=0

g(P_c)=M(P_c)-M_{\mathrm{target}}

Choose $P_c$ → integrate → locate the surface → return $M,R$ .

Adjust $P_c$ until $M(P_c)=M_{\rm target}$ .

Shooting wraps center-to-surface IVP integrations in an outer root solve.

Return to Project 4 and synthesize

Explain the mechanism

Halving $h$ halves endpoint error. Which order does this suggest?

An RK4 routine reuses $k_2$ as $k_3$ . Which state calculation did it skip?

Oscillator peaks arrive late despite small energy error. Which diagnostic detects this?

Discuss: Which checks apply to a static planet profile? Which are oscillator-specific?

State, energy, and phase measurements answer different questions.

Optional references

Euler’s finite-interval error bound

e_{n+1}\leq(1+hL)e_n+Ch^2

e_N\leq\frac{Ch}{L}(e^{LT}-1),\quad L>0

Exact initial data; fixed $T=Nh$ .
If $L=0$ , the bound is $CTh$ .

Smooth local defects plus Lipschitz propagation give first-order finite-time error. — Each computed state supplies the next slope; errors propagate.

Optional references

Why explicit midpoint is second order

y^{\prime\prime}=f_x+J_f f

$f_x=\partial f/\partial x$ : coordinate derivative.
$J_f=\partial f/\partial\mathbf y$ : state Jacobian.

k_2=f+\frac h2(f_x+J_f f)+O(h^2)

Multiply by $h$ ; add the original state.
The expansions agree through $h^2$ .

Stage-state and stage-coordinate shifts both contribute to the expansion.

Optional references

Implicit stability is not resolved accuracy

Backward Euler: first order. Global state error $O(h)$ .

u^\prime=-1000u,\qquad h=0.01

u_{n+1}=u_n+h f(t_{n+1},u_{n+1})

Amplification factor: $R=u_{n+1}/u_n$ .

R_{\rm Euler}=-9,\qquad R_{\rm BE}=1/11

Exact amplification: $e^{-10}$ .
Stable, but the fast transient is unresolved.

At the stated h=0.01, backward Euler stays positive but decays much more slowly than the exact fast transient. — Backward Euler stays stable while missing the rapid transient.

Optional references

Adaptive steps control an error estimate

Compare two approximations to estimate local error.
Accept or reject before advancing.
Internal steps need not match output spacing.

The error scale depends on component units and the requested observable.