This interactive tutorial compares numerical integration methods for solving ordinary differential equations (ODEs): Euler, Heun (Improved Euler), and Runge-Kutta 4th order (RK4). The test problem is a damped harmonic oscillator (spring-mass system with friction). You can choose linear stiffness (F = −kx), nonlinear stiffness (Cubic Stiff F = −kx³ or Cubic Soft Duffing-style), and damping type: linear (c·v), quadratic (c·v|v|, air resistance), or Coulomb (constant friction). The system is integrated as two first-order ODEs: dx/dt = v and dv/dt = restoring force + damping force. The simulation shows how step size (h) and algorithm choice affect accuracy and stability.

The visualization includes: (1) Spring view – up to three masses (Euler=red, Heun=orange, RK4=blue) on springs; (2) Phase portrait – velocity (v) vs position (x); (3) Time series – position x(t) vs time, with exact solution (green dashed) when linear and underdamped; (4) Global Error (vs Exact) – red = Euler error vs reference, blue = RK4 error vs reference (analytical for linear systems, numerical reference for nonlinear); (5) Step calculation – at the bottom, formula → numbers → result for Euler, Heun, and RK4 for the last step. A Numerical divergence box shows the Euclidean distance between Euler and RK4 in phase space (red when large).

Mathematical model

The second-order equation m x″ + (damping) + (restoring) = 0 is written as two first-order ODEs in the state vector y = [x, v]:

Restoring force (Nonlinear k): Linear −kx; Cubic Stiff −kx³; Cubic Soft −k(x − 0.1x³) (Duffing-style). Damping force (Damping Type): Linear −c v; Quadratic −c v|v|; Coulomb −c sign(v). The standard linear case is dv/dt = −(c/m)v − (k/m)x.

The three integrators advance the state y by step h as follows:

Euler is first-order (can gain energy and spiral outward when h is large), Heun is second-order, and RK4 is fourth-order — much more accurate and stable for oscillators.

Phase portrait: Plot (x, v). For damped systems the trajectory spirals inward to (0,0). Euler with large h can spiral outward; RK4 approximates the true spiral. Exact solution (linear, underdamped): ζ = c/(2√(km)) < 1, ω₀ = √(k/m), ω_d = ω₀√(1−ζ²); x(t) = e^−ζω₀t(A cos(ω_dt) + B sin(ω_dt)) with A = x₀, B = (v₀ + ζω₀x₀)/ω_d. For nonlinear systems, error is computed vs a numerical reference (high-precision RK4 micro-steps).

Simulation

The interactive simulator is below. Use the controls to explore the concepts described above.

 

Usage Example

Follow these steps to explore the Numerical Integration simulation:

1. Initial state: Load the page with default (damped) preset. You see: (1) Spring view – three masses (Euler=red, Heun=orange, RK4=blue) at the same initial position; (2) Phase portrait – velocity (v) vs position (x); (3) Time series – position x(t) vs time. Use Step Fwd or Run to advance; Step Bwd to go back; Reset to restart.
2. Compare methods: With Method: Euler + Heun + RK4, run the simulation. In the phase portrait, the damped system should spiral inward toward (0,0). Euler (red) may spiral outward when step size h is large – that is numerical instability (energy gain). RK4 (blue) stays close to the correct spiral.
3. Step size (h): Increase h (e.g. to 0.3–0.5). Euler diverges quickly; RK4 remains stable. Decrease h to see Euler “snap” closer to the curve. The Numerical divergence value (Euclidean distance between Euler and RK4 in phase space) grows when Euler drifts; it turns red when > 1.
4. Undamped (c = 0): Choose preset Undamped or set Damping (c) to 0. The exact solution is a circle in phase space. Euler spirals outward; RK4 stays near the circle. This shows that RK4 preserves energy much better for oscillators.
5. Euler unstable preset: Choose Euler unstable (large h) to see Euler explode quickly; the message “Lower the Step Size” appears. Click Reset and reduce h to continue.
6. Presets: Use the Preset dropdown: Default (damped), Undamped, Euler unstable, Stiff, Heavy damping. Each loads mass (m), damping (c), stiffness (k), step size (h), and initial x₀, v₀. Adjust sliders afterward; changing a slider does not change the preset name (you can reset and choose preset again to reapply).
7. Method: Switch to Euler only, Heun only, or RK4 only to focus on one method. The spring view and both graphs show only that method.

Parameters

• Step size (h): Time step for integration. Range: 0.01–0.5. Larger h makes Euler unstable; RK4 remains accurate for moderate h. All parameters apply live while Run is active.
• Mass (m), Damping (c), Stiffness (k): Physical parameters. ω₀ = √(k/m), ζ = c/(2√(km)). ζ < 1: underdamped; ζ = 1: critically damped; ζ > 1: overdamped.
• Damping Type: Linear (F_d = −c v), Quadratic (F_d = −c v|v|, air resistance), Coulomb (constant friction, zero at v=0).
• Nonlinear k: Linear (F = −kx), Cubic Stiff (F = −kx³), Cubic Soft (F = −k(x − 0.1x³), Duffing-style).
• x₀, v₀: Initial position and velocity at t = 0.
• Method: Euler only, Heun only, RK4 only, or Euler + Heun + RK4 (compare all).

Controls and visualizations

• Preset: Dropdown to load Default (damped), Undamped, Euler unstable (large h), Stiff, Heavy damping. Resets simulation when changed.
• Step Bwd / Step Fwd: Move one integration step backward or forward (stops Run if playing).
• Run / Stop: Start or stop continuous stepping.
• Reset: Restore initial conditions and clear history; use current parameter values.
• Numerical divergence: Euclidean distance between Euler and RK4 in (x, v) space; turns red when > 1.
• Spring view: Horizontal position of the mass(es) on a spring; red=Euler, orange=Heun, blue=RK4 when Method = all.
• Phase portrait: v vs x with grid; Auto Scale On/Off. Euler can spiral outward (unstable); RK4 spirals inward for damped systems.
• Time series: x(t) vs t with sliding time window. Green dashed = exact solution (linear, underdamped only). Red/orange/blue = Euler/Heun/RK4.
• Global Error (vs Exact): Red = Euler error vs reference, blue = RK4 error vs reference. For linear systems reference is analytical; for nonlinear (Cubic Stiff/Soft or Quadratic/Coulomb damping) reference is a high-precision numerical (RK4 micro-step) solution.
• Step calculation: Info box at the bottom shows for the last step: formula → numbers plugged in → result for Euler, Heun, and RK4 (including k₁…k₄ formulas for RK4).

Key concepts

• Damped harmonic oscillator: m x″ + damping + restoring = 0; state (x, v). Linear case: −(c/m)v − (k/m)x. Nonlinear stiffness (Cubic Stiff/Soft) and damping (Quadratic, Coulomb) expose algorithm differences clearly.
• Euler method: First-order; y_n+1 = y_n + h f(y_n). Tends to add energy in oscillators; diverges quickly for nonlinear or large h.
• Heun (Improved Euler): Second-order; k₁ = f(y_n), k₂ = f(y_n+h·k₁), y_n+1 = y_n + (h/2)(k₁+k₂). More accurate than Euler, less than RK4.
• RK4: Fourth-order; k₁…k₄ with mid/end samples; y_n+1 = y_n + (h/6)(k₁+2k₂+2k₃+k₄). Much better accuracy and stability for oscillators and nonlinear systems.
• Error reference: Linear (F=−kx, linear damping): error vs analytical exact solution. Nonlinear: error vs numerical reference (RK4 with 50 micro-steps per step) so the Global Error chart remains meaningful.
• Phase portrait: v vs x. Correct damped trajectory spirals inward; Euler with large h can spiral outward (numerical instability).
• Numerical divergence: Euclidean distance between Euler and RK4 in phase space; quantifies drift (fake energy).

Limitations

• One test problem. The integrators are compared on a single (damped harmonic / Duffing) oscillator. Behavior on stiff systems, chaotic ODEs, PDEs, or higher-dimensional state vectors can differ markedly.
• Three fixed-step methods. Only explicit Euler, Heun, and classical RK4 at a constant step h. Adaptive step-size control (RK45/Dormand–Prince), implicit/stiff solvers (backward Euler, BDF), and symplectic integrators are not included.
• No global error theory shown. The demo displays error curves but does not derive order-of-accuracy or stability-region (A-stability) analysis; "Lower the Step Size" is a heuristic, not a CFL/stability bound.
• Reference is itself numerical for nonlinear cases. Exact error is only available for the linear underdamped case; otherwise the "truth" is a fine-step RK4, so reported nonlinear errors are relative to an approximation.
• Floating-point and round-off. Very small h eventually trades truncation error for accumulated round-off; this regime is not explored.
• Scalar 1-DOF mechanics. A single mass on a spring in arbitrary units; multi-body coupling, forcing/resonance sweeps, and energy-conservation diagnostics beyond the divergence metric are out of scope.