Skip to content
Computational Lab.

RK4 vs Euler

The same oscillator, three integrators.

Intermediate / ODE / integration / error-analysis

The harmonic oscillator has an exact solution, so integration error is not a matter of opinion: forward Euler pumps energy into the system until it spirals away, semi-implicit Euler stays bounded but drifts in phase, and RK4 tracks the true curve to several decimal places. Change the timestep and watch the error orders move.

Tagged: RK4, Euler, symplectic, energy error, ODE

Live simulation

Lower dt to watch RK4 and semi-implicit Euler converge on the exact curve.

Speed

Lower dt to watch RK4 and semi-implicit Euler converge on the exact curve.

Mathematics

The harmonic oscillator θ̈ = −ω²θ is the standard test case for integrators because its exact solution, θ(t) = cos(ωt), is known. Every method's error is therefore measurable rather than estimated.

forward Euler: ω' = ω + a(θ)dt; θ' = θ + ω dt → energy grows, O(dt) globally

semi-implicit: ω' = ω + a(θ)dt; θ' = θ + ω' dt → symplectic, bounded energy

classical RK4: four stage evaluations, O(dt⁴) globally

Forward Euler is structurally wrong for oscillators: each step adds a little energy, so the amplitude spirals outward forever. Semi-implicit Euler fixes it by using the new velocity for the position update, which preserves a modified energy exactly and keeps the orbit bounded. RK4 simply converges four orders faster. The error metrics are computed from the energy itself, E = ½ω²θ² + ½ω̇², relative to E₀.

How it works

All three integrators advance the same initial condition in lockstep, and the plot draws the exact solution as a reference line beneath them. The history buffer keeps a rolling window so the growing amplitude of Euler's spiral is visible without the axes ever moving.

Cut dt from 0.02 to 0.005 and watch RK4's error drop by roughly 256×, which is what fourth-order convergence predicts: 4⁴ = 256.

Implementation

src/engine/experiments/rk4-vs-euler.ts · 167 lines · the module running above

TypeScript · fixed-timestep engine

1import { BaseExperiment } from "../core/base-experiment";
2import type { ParameterDef, ThemeColors } from "../core/types";
3
4const SAMPLES = 900;
5
6/**
7 * One oscillator, three integrators. theta'' = -omega^2 * theta has the exact
8 * solution theta(t) = cos(omega t), so every method's energy error is a real
9 * number you can read off the plot: forward Euler spirals out, semi-implicit
10 * Euler stays bounded with phase error, RK4 stays flat.
11 */
12export class IntegratorComparisonExperiment extends BaseExperiment {
13 readonly id = "rk4-vs-euler";
14
15 private dt = 0.02;
16 private omega = 1;
17 private time = 0;
18 private steps = 0;
19 private euler = { theta: 1, omegaDot: 0 };
20 private semi = { theta: 1, omegaDot: 0 };
21 private rk4 = { theta: 1, omegaDot: 0 };
22 private history: { euler: number; semi: number; rk4: number; exact: number }[] = [];
23 private errorEuler = 0;
24 private errorSemi = 0;
25 private errorRk4 = 0;
26
27 protected params(): ParameterDef[] {
28 return [
29 { key: "dt", label: "Timestep dt", min: 0.005, max: 0.08, step: 0.005, defaultValue: 0.02, unit: "s" },
30 { key: "omega", label: "Frequency ω", min: 0.5, max: 3, step: 0.1, defaultValue: 1, unit: "rad/s" },
31 { key: "steps", label: "Steps per tick", min: 1, max: 40, step: 1, defaultValue: 8 },
32 ];
33 }
34
35 protected onReset(): void {
36 this.euler = { theta: 1, omegaDot: 0 };
37 this.semi = { theta: 1, omegaDot: 0 };
38 this.rk4 = { theta: 1, omegaDot: 0 };
39 this.history = [];
40 this.time = 0;
41 this.steps = 0;
42 this.errorEuler = 0;
43 this.errorSemi = 0;
44 this.errorRk4 = 0;
45 }
46
47 protected onParameterChange(key: string): void {
48 if (key !== "steps") this.reset();
49 }
50
51 /** Energy relative error |E/E0 - 1| for theta'' = -w^2 theta. */
52 private relativeError(theta: number, omegaDot: number): number {
53 const energy = 0.5 * (omegaDot * omegaDot + this.omega * this.omega * theta * theta);
54 const e0 = 0.5 * this.omega * this.omega;
55 return Math.abs(energy / e0 - 1);
56 }
57
58 private accel(theta: number): number {
59 return -this.omega * this.omega * theta;
60 }
61
62 protected onUpdate(): void {
63 const dt = this.num("dt");
64 this.omega = this.num("omega");
65 const stepsPerTick = this.num("steps");
66 for (let s = 0; s < stepsPerTick; s++) {
67 // Forward Euler: unstable for oscillators (energy grows).
68 const eulerA = this.accel(this.euler.theta);
69 this.euler.theta += this.euler.omegaDot * dt;
70 this.euler.omegaDot += eulerA * dt;
71 // Semi-implicit (symplectic) Euler: bounded, phase-shifted.
72 this.semi.omegaDot += this.accel(this.semi.theta) * dt;
73 this.semi.theta += this.semi.omegaDot * dt;
74 // Classical RK4.
75 const a1 = this.accel(this.rk4.theta);
76 const b1 = this.accel(this.rk4.theta + (this.rk4.omegaDot * dt) / 2);
77 const c1 = this.accel(this.rk4.theta + (this.rk4.omegaDot * dt) / 2 + (a1 * dt * dt) / 4);
78 const d1 = this.accel(this.rk4.theta + this.rk4.omegaDot * dt + (b1 * dt * dt) / 2);
79 const nw =
80 this.rk4.omegaDot + (dt / 6) * (a1 + 2 * b1 + 2 * c1 + d1);
81 const nt =
82 this.rk4.theta +
83 (dt / 6) *
84 (this.rk4.omegaDot +
85 2 * (this.rk4.omegaDot + a1 * dt / 2) +
86 2 * (this.rk4.omegaDot + b1 * dt / 2) +
87 (this.rk4.omegaDot + c1 * dt));
88 this.rk4.omegaDot = nw;
89 this.rk4.theta = nt;
90
91 this.time += dt;
92 this.steps += 1;
93 this.history.push({
94 euler: this.euler.theta,
95 semi: this.semi.theta,
96 rk4: this.rk4.theta,
97 exact: Math.cos(this.omega * this.time),
98 });
99 if (this.history.length > SAMPLES) this.history.shift();
100 }
101 this.errorEuler = this.relativeError(this.euler.theta, this.euler.omegaDot);
102 this.errorSemi = this.relativeError(this.semi.theta, this.semi.omegaDot);
103 this.errorRk4 = this.relativeError(this.rk4.theta, this.rk4.omegaDot);
104 }
105
106 render(ctx: CanvasRenderingContext2D, theme: ThemeColors): void {
107 ctx.fillStyle = theme.bg;
108 ctx.fillRect(0, 0, this.width, this.height);
109 const pad = 18;
110 const mid = this.height / 2;
111 const scaleY = (this.height / 2 - pad) / 2.2;
112 // Zero line
113 ctx.strokeStyle = theme.line;
114 ctx.lineWidth = 1;
115 ctx.beginPath();
116 ctx.moveTo(0, mid);
117 ctx.lineTo(this.width, mid);
118 ctx.stroke();
119 const plot = (pick: (h: (typeof this.history)[number]) => number, color: string, width = 1.5) => {
120 if (this.history.length < 2) return;
121 ctx.strokeStyle = color;
122 ctx.lineWidth = width;
123 ctx.beginPath();
124 for (let i = 0; i < this.history.length; i++) {
125 const x = (i / (SAMPLES - 1)) * this.width;
126 const y = mid - pick(this.history[i]) * scaleY;
127 if (i === 0) ctx.moveTo(x, y);
128 else ctx.lineTo(x, y);
129 }
130 ctx.stroke();
131 };
132 plot((h) => h.exact, theme.fgTertiary, 1);
133 plot((h) => h.euler, theme.viz[3]);
134 plot((h) => h.semi, theme.viz[1]);
135 plot((h) => h.rk4, theme.accent, 2);
136 // Legend
137 ctx.font = "11px ui-monospace, monospace";
138 ctx.fillStyle = theme.fgSecondary;
139 ctx.fillText("EXACT", 10, 16);
140 ctx.fillStyle = theme.viz[3];
141 ctx.fillText("EULER", 62, 16);
142 ctx.fillStyle = theme.viz[1];
143 ctx.fillText("SEMI-IMPLICIT", 112, 16);
144 ctx.fillStyle = theme.accent;
145 ctx.fillText("RK4", 216, 16);
146 }
147
148 getMetrics() {
149 const fmt = (e: number) => e.toExponential(2);
150 return {
151 simTime: `${this.time.toFixed(2)} s`,
152 steps: this.steps,
153 eulerError: fmt(this.errorEuler),
154 semiError: fmt(this.errorSemi),
155 rk4Error: fmt(this.errorRk4),
156 };
157 }
158
159 describe(): string {
160 return `After ${this.time.toFixed(1)} simulated seconds, RK4's energy error is ${this.errorRk4.toExponential(1)} while forward Euler has drifted to ${this.errorEuler.toExponential(1)}.`;
161 }
162
163 entities(): number {
164 return this.history.length * 4;
165 }
166}
167

Performance

Time per tick
O(steps per tick); RK4 costs four evaluations per step
Space
A rolling 900-sample history per curve
Notes
This page is the numerical-methods category's thesis: the integrator is part of the physics whether you choose it or not.

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.