Runge-Kutta Simulator Back
Numerical Methods Simulator

Runge-Kutta Simulator — RK4 vs Euler Accuracy Comparison

Integrate the ODE dy/dx = -k*y in real time using both classical fourth-order Runge-Kutta (RK4) and the explicit Euler method, then compare against the exact decay y = y0*exp(-k*x). Final-value statistics and Euler relative error appear instantly, and the cumulative absolute error is plotted on a log10 axis so that the concept of convergence order becomes tangible.

Parameters
Decay rate k
1/x
Step size h
End point x_end
Initial value y0

Defaults: k=0.5, h=0.1, x_end=4, y0=10. The exact value is y(4) = 10*e^-2 ≈ 1.3534. Increasing h explodes the Euler error.

Results
RK4 final y(x_end)
Euler final y(x_end)
Exact y(x_end)
Euler relative error
Solution comparison (Exact / RK4 / Euler)

x axis: x in [0, x_end]; y axis: y. White heavy line: exact y = y0*exp(-k*x); blue: RK4; red: Euler. With small h, RK4 nearly overlaps the exact curve while Euler droops below it; with large h, the Euler gap becomes obvious.

Cumulative error |y_num - y_exact| (log10 scale)

x axis: x; y axis: log10|y_num - y_exact|. Blue: RK4 absolute error; red: Euler absolute error. RK4 is several decades smaller than Euler (about 1e5 difference at h=0.1), and both grow with x at very different rates.

Theory & Key Formulas

The ODE and its exact solution:

$$\frac{dy}{dx} = -k\,y,\qquad y(0)=y_0 \;\Longrightarrow\; y(x) = y_0\,e^{-kx}$$

Explicit Euler (first-order, global error O(h)):

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

Classical fourth-order Runge-Kutta (RK4, global error O(h^4)):

$$\begin{aligned}k_1 &= h\,f(x_n, y_n)\\ k_2 &= h\,f(x_n+\tfrac{h}{2}, y_n+\tfrac{k_1}{2})\\ k_3 &= h\,f(x_n+\tfrac{h}{2}, y_n+\tfrac{k_2}{2})\\ k_4 &= h\,f(x_n+h, y_n+k_3)\\ y_{n+1} &= y_n + \tfrac{1}{6}(k_1+2k_2+2k_3+k_4)\end{aligned}$$

$k$ is the decay rate (time constant $\tau=1/k$), $h$ is the step size and $y_0$ is the initial value. Euler uses only the starting slope, RK4 takes a Simpson-style weighted average of four slopes; halving $h$ shrinks Euler error by 1/2 and RK4 error by 1/16.

What is the Runge-Kutta Simulator

🙋
My textbook keeps saying RK4 is more accurate than Euler, but how much more accurate, really? I would love to see the numbers.
🎓
Try the defaults of this tool: k=0.5, h=0.1, x_end=4, y0=10. Solving dy/dx = -0.5*y to x=4 gives an exact value of 10*e^-2 = 1.3534. The result cards show Euler at 1.2851 (about -5.05% off) and RK4 at 1.3534, matching the exact value to four decimals. On the log-scale error plot below, RK4 error sits around 1e-5 while Euler error is around 1e-1, a factor of ten thousand difference.
🙋
Ten thousand times more accurate, just like that? Both methods are integrating the same ODE.
🎓
Both yes, but they evaluate the slope a different number of times per step. Euler uses only the slope at the start of the interval, so anything past the first order in the Taylor expansion becomes error: local error O(h^2), global error O(h). RK4 evaluates the slope at the start, twice at the midpoint, and at the end, then averages with Simpson-style weights (1, 2, 2, 1)/6. That matches the Taylor expansion through fourth order, giving local error O(h^5) and global error O(h^4). At h=0.1 that is 1e-4 instead of 1e-1.
🙋
When I pressed play to sweep h, the Euler curve started oscillating around h=0.5 and finally blew up. RK4 looked fine.
🎓
That is numerical stability. Euler is only stable when |1 - kh| <= 1, so for k=0.5 the limit is h=4, but oscillation appears well before that. RK4 has a much larger stability region, roughly |λh| <= 2.78. On stiff problems Euler needs ridiculously small steps to stay stable, which is why production CAE solvers reach for implicit schemes such as Hilber-Hughes-Taylor (Abaqus/Standard) or central differences (LS-DYNA) depending on the application.
🙋
So RK4 is always the winner? Should I just default to it?
🎓
Not quite. RK4 evaluates f four times per step, so its cost per step is four times Euler. For an apples-to-apples comparison you would let Euler use h/4 against RK4 with h, and RK4 still wins by orders of magnitude. In practice the standard is Dormand-Prince RK45 with built-in step-size control: this is what MATLAB ode45, SciPy RK45 and Mathematica NDSolve all use. RK4 is the classical textbook reference and the gateway to understanding all of them.

FAQ

The Euler method estimates the slope across the interval [x_n, x_n+h] using only the value at the start, so its local error is O(h^2) and its global error is just O(h). RK4 evaluates the slope at four points (start, midpoint twice, end) and averages them with Simpson-like weights (1, 2, 2, 1)/6, which gives local error O(h^5) and global error O(h^4). With the defaults of this tool (k=0.5, h=0.1, x_end=4, y0=10) Euler returns y(4)=1.2851, which is 5.05% below the exact 1.3534, while RK4 returns 1.3534 with an error well below 1e-9. To match RK4 accuracy with Euler you would need to shrink h by a factor of ten thousand, raising the cost thousands of times.
Larger h means fewer steps and faster runs, but the local error scales as a power of h. Euler error grows linearly with h, while RK4 error grows as h^4, so doubling h doubles the Euler error but multiplies the RK4 error by sixteen. For stiff problems such as dy/dx = -k*y, the Euler solution diverges once |1 - kh| > 1 (for k=0.5 that means h > 4). Pressing play to sweep h from 0.005 up to 1.0 lets you watch Euler start oscillating near h=0.5 and break down completely beyond h=2, while RK4 remains stable up to about h = 2.78, which is why RK4 is still usable on mildly stiff problems.
Common ODE solvers include (1) explicit Euler (first order, simplest), (2) the midpoint rule (second order), (3) Heun's or improved Euler method (second order), (4) RK4 (fourth order, used here), (5) Adams-Bashforth multistep schemes that reuse past values to save evaluations, (6) backward Euler, the trapezoidal rule, and BDF for stiff problems, and (7) Dormand-Prince RK45, the adaptive scheme used by MATLAB ode45 and SciPy RK45. In CAE practice LS-DYNA uses central differences for explicit dynamics and Abaqus/Standard uses Hilber-Hughes-Taylor (the alpha method) for implicit dynamics. RK4 is the textbook classic, but production codes generally rely on adaptive RK45 with built-in error control.
This tool integrates only the scalar linear ODE dy/dx = -k*y and ignores (1) coupled ODE systems such as Lorenz attractors, (2) nonlinear ODEs such as Van der Pol or Lotka-Volterra, (3) time integration of partial differential equations such as heat or wave equations, (4) adaptive step-size control like RK45, (5) automatic stiffness detection and switching to implicit solvers, and (6) symplectic integration that preserves energy in Hamiltonian systems. Real CAE solvers combine these ideas through schemes such as Newmark-beta for structural dynamics, Crank-Nicolson for CFD time stepping, and Verlet integration for orbital mechanics. This tool is best for building intuition about convergence order.

Real-world applications

Structural dynamics (FEM transient analysis): seismic and crash analyses time-integrate M*u'' + C*u' + K*u = F(t). Abaqus/Standard uses the Hilber-Hughes-Taylor alpha method (implicit, second-order accurate), while LS-DYNA uses the central difference method (explicit, second order). Explicit schemes avoid stiffness matrix inversion and excel at crash simulations, while implicit ones tolerate much larger time steps for long-duration response. That trade-off between accuracy, stability and cost is exactly what you feel by sliding h in this tool. Newmark-beta with beta=1/4 and gamma=1/2 is unconditionally stable and is the most widely used integrator in commercial CAE.

CFD time stepping: OpenFOAM icoFoam uses implicit Euler or Crank-Nicolson, while SU2 ships RK4 as a standard option. The Navier-Stokes equations are inherently nonlinear and stiff, so the CFL (Courant-Friedrichs-Lewy) condition limits the time step. Watching Euler diverge as h grows in this tool is the same phenomenon as a CFD solver blowing up when CFL exceeds one. Large-eddy simulations (LES) typically reach for TVD-RK3 (third-order Runge-Kutta with Total Variation Diminishing) to suppress oscillations around shocks.

Orbital mechanics (satellites and planetary motion): NASA's SPICE Toolkit and GMAT use RK7(8) or the Bulirsch-Stoer method for satellite orbits. Low-Earth orbit work demands accuracy on the order of 0.1 second per revolution, and decade-long propagations rely on Verlet integration or symplectic schemes to preserve energy. The exponential decay y = y0*exp(-k*x) handled here has a closed form, just as the Kepler two-body problem does, and the short-term vs long-term trade-off between RK4 and symplectic integrators mirrors the same intuition.

Control engineering (real-time integration): the integral term of a PID controller, the simulation of a state-space model x' = Ax + Bu, and the prediction step of a Kalman filter are all ODE-integration tasks. Real-time controllers often run fixed Euler at a 1 ms step (cheap on FPGAs), while offline simulators such as CARLA or MATLAB Simulink default to ode45 (RK45). Sliding k and h here, you can verify the rule of thumb that the step size should be at least ten times finer than the control bandwidth (k*h much less than 0.1) for safe operation.

Common misconceptions and caveats

The most common misconception is that RK4 is always more accurate than anything else. The global error of RK4 is bounded by C*h^4, but the constant C depends on the higher derivatives of the solution. For highly oscillatory or non-smooth right-hand sides (such as those involving Heaviside steps), C grows large and the apparent accuracy collapses. Crash analyses, for instance, prefer the LeapFrog scheme or TVD-RK methods that preserve wave fronts. RK4 is excellent on smooth exponential decays like the one in this tool, but always match the method to the structure of your problem.

Another popular myth is that smaller step sizes are always better. Shrinking h reduces discretization error, but it inflates the total number of operations and accumulates floating-point round-off. Double-precision (IEEE 754) has a machine epsilon of about 2.2e-16, so once h drops below 1e-8 round-off begins to dominate and accuracy actually degrades (the optimal h is near sqrt(2*epsilon), about 1e-8). The minimum h of this tool, 0.005, stays well above that limit; edit the source if you want to probe smaller values.

Finally, do not assume that RK4 handles stiff problems. The RK4 stability region is bounded roughly by |lambda*h| <= 2.78. A chemical system with two species at time constants of 1 s and 1 us would require h far below 1 us, making the computation infeasible at the second-time scale of interest. Stiff problems demand implicit integrators (backward Euler, BDF, Rosenbrock), embodied by MATLAB ode15s and SciPy BDF and Radau. That is why every CAE solver lets you choose between explicit and implicit modes.