Skip to content
Computational Lab.

Adaptive Step Size (RK45)

Small steps where it matters, big steps where it doesn't.

Advanced / ODE / integration / adaptive

The Dormand-Prince 5(4) embedded pair computes a fifth-order step and a fourth-order estimate from the same stages; their difference is an error estimate, and the controller shrinks or grows the step to hold that estimate near tolerance. Watch the step trace collapse through the stiff transient and relax afterward - then compare the accepted step count against fixed-step RK4 at dt = 1/60.

Tagged: Dormand-Prince, embedded pair, step-size control, stiff transient

Live simulation

Each mark in the lower band is one accepted step; watch h collapse in the transient.

Speed

Each mark in the lower band is one accepted step; watch h collapse in the transient.

Mathematics

An embedded Runge-Kutta pair computes two solutions of different order from the same stage evaluations. Dormand-Prince's 5(4) pair gives seventh stages producing a fifth-order step and a fourth-order estimate; their difference estimates the local error, and a controller picks the next step size to hold that estimate near tolerance.

err = |dt · Σ (b_i − b̂_i) k_i|

controller: h ← h · clamp(0.9 (tol / err)^{1/5}, 0.2, 5)

test problem: y' = λ (sin t − y) + cos t, exact: y(t) = sin t

The test problem's transient e^{−λt} is where the controller earns its keep: with λ = 30, the first fraction of a second demands steps orders of magnitude smaller than the smooth sin tracking that follows. The step trace shows exactly that collapse and recovery. Because the exact solution is known, the reported max error is the real error, and the fixed-step comparison (dt = 1/60, 480 steps per pass) shows what adaptivity buys.

How it works

Each tick performs a fixed number of accepted steps. A step is attempted, rejected and retried with a shrunk dt while its embedded error exceeds tolerance, then accepted with a controller-suggested next h. The trace stores t, y, and h per accepted step; one full pass over the 8-second horizon snapshots the metrics before restarting.

Everything is deterministic: same parameters, same step sequence - the controller is part of the simulation, not hidden machinery.

Implementation

src/engine/experiments/adaptive-rk45.ts · 226 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import type { ParameterDef, ThemeColors } from "../core/types";
3
4const HORIZON = 8;
5
6/**
7 * Adaptive step-size integration with the Dormand-Prince 5(4) embedded pair.
8 * The ODE y' = lambda (sin t - y) + cos t has exact solution y = sin t, so
9 * the controller's tolerance is checked against the truth, and the step
10 * trace shows h shrinking through the stiff transient and relaxing after.
11 */
12export class AdaptiveRk45Experiment extends BaseExperiment {
13 readonly id = "adaptive-rk45";
14
15 private t = 0;
16 private y = 0;
17 private h = 0.05;
18 private lambda = 30;
19 private accepted = 0;
20 private rejected = 0;
21 private maxError = 0;
22 private lastPass = { accepted: 0, rejected: 0, maxError: 0 };
23 private trace: { t: number; y: number; h: number }[] = [];
24
25 protected params(): ParameterDef[] {
26 return [
27 { key: "lambda", label: "Stiffness lambda", min: 5, max: 60, step: 5, defaultValue: 30 },
28 { key: "tolExp", label: "Tolerance 10^", min: -10, max: -3, step: 1, defaultValue: -6 },
29 { key: "steps", label: "Steps per tick", min: 1, max: 40, step: 1, defaultValue: 8 },
30 ];
31 }
32
33 protected onReset(): void {
34 this.t = 0;
35 this.y = 0;
36 this.h = 0.05;
37 this.accepted = 0;
38 this.rejected = 0;
39 this.maxError = 0;
40 this.trace = [];
41 void this.lastPass;
42 this.lastPass = { accepted: 0, rejected: 0, maxError: 0 };
43 }
44
45 protected onParameterChange(key: string): void {
46 if (key !== "steps") this.reset();
47 }
48
49 private derivative(t: number, y: number): number {
50 return this.lambda * (Math.sin(t) - y) + Math.cos(t);
51 }
52
53 /** One Dormand-Prince attempt of size dt from (t, y): the 5th-order
54 * solution increment and the embedded 4th-order error estimate. */
55 private dpTry(dt: number): { advance: number; err: number } {
56 const k1 = this.derivative(this.t, this.y);
57 const k2 = this.derivative(this.t + dt * (1 / 5), this.y + dt * (1 / 5) * k1);
58 const k3 = this.derivative(
59 this.t + dt * (3 / 10),
60 this.y + dt * ((3 / 40) * k1 + (9 / 40) * k2),
61 );
62 const k4 = this.derivative(
63 this.t + dt * (4 / 5),
64 this.y + dt * ((44 / 45) * k1 - (56 / 15) * k2 + (32 / 9) * k3),
65 );
66 const k5 = this.derivative(
67 this.t + dt * (8 / 9),
68 this.y +
69 dt * ((19372 / 6561) * k1 - (25360 / 2187) * k2 + (64448 / 6561) * k3 - (212 / 729) * k4),
70 );
71 const k6 = this.derivative(
72 this.t + dt,
73 this.y +
74 dt *
75 ((9017 / 3168) * k1 -
76 (355 / 33) * k2 +
77 (46732 / 5247) * k3 +
78 (49 / 176) * k4 -
79 (5103 / 18656) * k5),
80 );
81 const advance =
82 dt *
83 ((35 / 384) * k1 +
84 (500 / 1113) * k3 +
85 (125 / 192) * k4 -
86 (2187 / 6784) * k5 +
87 (11 / 84) * k6);
88 const k7 = this.derivative(this.t + dt, this.y + advance);
89 const errCoefficient =
90 (5179 / 57600) * k1 +
91 (7571 / 16695) * k3 +
92 (393 / 640) * k4 -
93 (92097 / 339200) * k5 +
94 (187 / 2100) * k6 +
95 (1 / 40) * k7 -
96 ((35 / 384) * k1 +
97 (500 / 1113) * k3 +
98 (125 / 192) * k4 -
99 (2187 / 6784) * k5 +
100 (11 / 84) * k6);
101 return { advance, err: Math.abs(dt * errCoefficient) };
102 }
103
104 protected onUpdate(): void {
105 const tolExp = this.num("tolExp");
106 const tol = Math.pow(10, tolExp);
107 const stepsBudget = this.num("steps");
108 this.lambda = this.num("lambda");
109 for (let s = 0; s < stepsBudget && this.t < HORIZON; s++) {
110 let dt = this.h;
111 let attempt = this.dpTry(dt);
112 let attempts = 0;
113 while (attempt.err > tol && attempts < 30 && dt > 1e-7) {
114 this.rejected += 1;
115 dt *= Math.max(0.2, 0.9 * Math.pow(tol / attempt.err, 1 / 5));
116 attempt = this.dpTry(dt);
117 attempts += 1;
118 }
119 this.y += attempt.advance;
120 this.t += dt;
121 this.accepted += 1;
122 this.maxError = Math.max(this.maxError, Math.abs(this.y - Math.sin(this.t)));
123 this.trace.push({ t: this.t, y: this.y, h: dt });
124 if (this.trace.length > 3000) this.trace.shift();
125 // Standard step-size controller for the next attempt.
126 this.h = Math.min(
127 0.5,
128 Math.max(
129 1e-6,
130 dt * Math.min(5, Math.max(0.2, 0.9 * Math.pow(tol / Math.max(attempt.err, 1e-14), 1 / 5))),
131 ),
132 );
133 }
134 if (this.t >= HORIZON) {
135 this.lastPass = { accepted: this.accepted, rejected: this.rejected, maxError: this.maxError };
136 this.t = 0;
137 this.y = 0;
138 this.h = 0.05;
139 this.accepted = 0;
140 this.rejected = 0;
141 this.maxError = 0;
142 this.trace = [];
143 }
144 }
145
146 render(ctx: CanvasRenderingContext2D, theme: ThemeColors): void {
147 ctx.fillStyle = theme.bg;
148 ctx.fillRect(0, 0, this.width, this.height);
149 const pad = 16;
150 const solutionBottom = this.height - 84;
151 const toX = (t: number) => pad + (t / HORIZON) * (this.width - 2 * pad);
152 const yScale = (this.height - 140) / 2.6;
153 const toY = (y: number) => solutionBottom / 2 - y * yScale + solutionBottom * 0;
154 const midY = solutionBottom / 2;
155 ctx.strokeStyle = theme.line;
156 ctx.beginPath();
157 ctx.moveTo(pad, midY);
158 ctx.lineTo(this.width - pad, midY);
159 ctx.stroke();
160 // Exact solution
161 ctx.strokeStyle = theme.fgTertiary;
162 ctx.setLineDash([4, 4]);
163 ctx.beginPath();
164 for (let px = 0; px <= this.width - 2 * pad; px += 2) {
165 const t = (px / (this.width - 2 * pad)) * HORIZON;
166 const y = midY - Math.sin(t) * yScale;
167 if (px === 0) ctx.moveTo(pad + px, y);
168 else ctx.lineTo(pad + px, y);
169 }
170 ctx.stroke();
171 ctx.setLineDash([]);
172 // Adaptive solution + step trace
173 if (this.trace.length > 1) {
174 ctx.strokeStyle = theme.accent;
175 ctx.lineWidth = 2;
176 ctx.beginPath();
177 for (let i = 0; i < this.trace.length; i++) {
178 const point = this.trace[i];
179 if (i === 0) ctx.moveTo(toX(point.t), midY - point.y * yScale);
180 else ctx.lineTo(toX(point.t), midY - point.y * yScale);
181 }
182 ctx.stroke();
183 // Step-size band: one tick per accepted step, height proportional to h.
184 const bandBase = this.height - 44;
185 ctx.strokeStyle = theme.viz[1];
186 ctx.lineWidth = 1;
187 ctx.beginPath();
188 for (const point of this.trace) {
189 const x = toX(point.t);
190 const h = Math.min(34, (Math.log10(point.h / 1e-5) / Math.log10(0.5 / 1e-5)) * 34);
191 ctx.moveTo(x, bandBase);
192 ctx.lineTo(x, bandBase - h);
193 }
194 ctx.stroke();
195 ctx.fillStyle = theme.fgTertiary;
196 ctx.font = "11px ui-monospace, monospace";
197 ctx.fillText("STEP SIZE h(t)", pad, this.height - 30);
198 ctx.strokeStyle = theme.line;
199 ctx.beginPath();
200 ctx.moveTo(pad, bandBase);
201 ctx.lineTo(this.width - pad, bandBase);
202 ctx.stroke();
203 }
204 void toY;
205 }
206
207 getMetrics() {
208 return {
209 lambda: this.lambda,
210 tolerance: `1e${this.num("tolExp")}`,
211 adaptiveSteps: this.lastPass.accepted,
212 rejected: this.lastPass.rejected,
213 fixedRk4Steps: HORIZON * 60,
214 maxError: this.lastPass.maxError.toExponential(2),
215 };
216 }
217
218 describe(): string {
219 return `Adaptive Dormand-Prince covered the stiff transient in ${this.lastPass.accepted} accepted steps where fixed dt=1/60 RK4 needs ${HORIZON * 60}; max error ${this.lastPass.maxError.toExponential(1)}.`;
220 }
221
222 entities(): number {
223 return this.trace.length;
224 }
225}
226

Performance

Time per tick
7 derivative evaluations per attempted step
Space
The t/y/h trace, 3000 samples
Notes
Tighten the tolerance slider and the step count rises; loosen it and the steps grow - the controller responds exactly as the theory says.

Live numbers (FPS, frame time, steps per second, microseconds per step, entity count) are in the performance overlay on the simulation above: the pulse icon in its top-right corner.