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).

NOTE: Use Step Bwd / Step Fwd to move one step backward or forward; Run/Stop to animate. If the simulation “explodes” (values go to infinity), the message Lower the Step Size appears – reduce h and click Reset. With Damping (c) = 0 (undamped), the exact solution is a perfect circle in phase space; Euler tends to spiral outward while RK4 stays near the circle.

Mathematical model

Second-order form: m x″ + (damping) + (restoring) = 0. Let v = x′. Then dx/dt = v and dv/dt = (restoring force + damping force)/m.

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) (constant magnitude, zero at v=0).

Standard linear case: dv/dt = −(c/m)v − (k/m)x with mass m, damping c, stiffness k.

Numerical methods (state vector y = [x, v]):

• Euler: y_n+1 = y_n + h·f(t_n, y_n). First-order; can gain energy and spiral outward when h is large.
• Heun (Improved Euler): k₁ = f(y_n), k₂ = f(y_n + h·k₁), y_n+1 = y_n + (h/2)(k₁ + k₂). Second-order.
• RK4: Four slope samples k₁,…,k₄; y_n+1 = y_n + (h/6)(k₁ + 2k₂ + 2k₃ + k₄). 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).

 

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.

Tip: To see Euler’s “death spiral” clearly: set Damping (c) to 0, Step size (h) to about 0.3, then Run. The red (Euler) trajectory spirals outward while the blue (RK4) trajectory stays near a circle. The Numerical divergence value increases as Euler gains fake energy.

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).