Visualize the absolute stability of the numerical integrators used to solve ordinary differential equations (Euler, RK2, RK4). Change the step size h and the eigenvalue λ to see the amplification factor R(z), the complex-plane stability region and whether the numerical solution diverges, all in real time.
Parameters
Integration method
Sets the order of the stability function R(z)
Test equation λ
Eigenvalue of y'=λy. λ<0 means the true solution decays
Step size h
Time interval advanced per step
Integration time t_end
s
Total time over which to integrate
Results
—
z = hλ
—
|R(z)| amplification
—
Stability verdict
—
Numerical final value
—
Exact final value
—
Maximum error
—
Complex-plane absolute stability region |R(z)| ≤ 1
The horizontal axis is the real part of z, the vertical axis the imaginary part. The shaded area is the stability region of the selected method. The current operating point z=hλ (on the real axis) is green = stable / red = unstable.
The test equation y'=λy (exact solution y(t)=y₀e^(λt)). One step advances by simply multiplying by the stability function R(z), so after n steps y_n = R(z)ⁿ·y₀.
The stability functions of each method. RK2 is $R=1+z+z^{2}/2$. The higher the order, the closer R(z) approaches the exponential eᶻ and the larger the stability region.
The condition for absolute stability. The set of z that satisfies it is the "absolute stability region". For a decaying system (λ<0) choose h so that z=hλ falls inside it.
What is Runge-Kutta Stability?
🙋
I learned that when you solve a differential equation on a computer, the smaller the step size h, the more accurate it gets. So if I just keep h small, I'm safe, right?
🎓
That is true as far as accuracy goes. But "accuracy" and "stability" are two different things. Stability is about the phenomenon where, if h is too large, the numerical solution diverges further and further and flies off in a completely wrong direction. For example, the true solution may quietly decay toward zero, yet the numerical solution alone oscillates and blows up to infinity. Making h smaller does fix it, but if you do not know "how large h can be", your computation just gets needlessly heavy.
🙋
What decides whether it diverges? When I dragged λ very negative on the left, it suddenly turned red and said "unstable".
🎓
Good catch. The key is a single dimensionless quantity, z = hλ. Solve one step of the test equation y'=λy and the numerical solution always takes the form y_{n+1} = R(z)·y_n. This R(z) is called the "stability function". Advance n steps and you get y_n = R(z)ⁿ. So if |R(z)| is larger than 1, the value swells with every step and diverges. If |R(z)| ≤ 1, it decays. Pushing λ negative moved z=hλ outside the stability region — that is why it turned red.
🙋
I see! So the difference between the stability regions of Euler and RK4 is what I'm seeing in that complex-plane plot above?
🎓
Exactly. Forward Euler has R(z)=1+z, so its stability region is a small unit disk centred at -1. On the real axis only z∈[-2,0] is stable. RK4 has R(z)=1+z+z²/2+z³/6+z⁴/24, giving a much larger leaf-shaped region that reaches about z≈-2.78 on the real axis. So for the same λ, RK4 lets you take a larger step size h. Choosing the method wisely is itself a way to cut computational cost.
🙋
Then surely it always pays to use RK4 and push h to the maximum?
🎓
I understand the temptation, but that is the trap. Right at the edge of the stability region you indeed do not diverge, but the accuracy collapses. When |R(z)| is close to 1, a solution that should decay fast lingers forever. Stability is the minimum condition for not blowing up; accuracy is how close to the truth you are. In practice engineers often use about half the stability-limit h to keep an accuracy margin. Especially for stiff equations — systems that mix fast and slow components — this constraint is harsh for explicit methods, and switching to an implicit method like backward Euler is the standard move.
🙋
Where do stiff equations actually show up in practice?
🎓
Chemical-reaction simulation is the classic example. Fast reactions reach equilibrium in milliseconds, while slow reactions take hours. Try to solve both in one system and the fast component's stability condition forces you to shrink h drastically. Circuit analysis, combustion, and drug pharmacokinetics in the body are the same. Drag λ to -30 in this tool and you will see that an explicit method needs h down around 0.07 to stay stable. That is the severity of stiffness.
Frequently Asked Questions
Absolute stability is the property that, when the test equation y'=λy (with λ<0, so the true solution decays) is solved numerically, the numerical solution also keeps decaying instead of blowing up. If R(z) is the per-step amplification factor (the stability function, z=hλ), then after n steps the numerical solution is y_n = R(z)^n. Only when |R(z)| ≤ 1 does the solution stay bounded as steps accumulate, which is absolute stability. When |R(z)| > 1 the numerical solution diverges exponentially even though the true solution is decaying.
On the real axis, forward Euler R(z)=1+z is stable only for z∈[-2,0] (a unit disk centred at -1 in the complex plane). RK4 R(z)=1+z+z²/2+z³/6+z⁴/24 is stable up to about z∈[-2.78,0], and its stability region is a larger lobe. So for the same decay rate λ, RK4 tolerates a larger step size h. Where Euler diverges once h exceeds roughly 2/|λ|, RK4 keeps going until about 2.78/|λ|.
For stability, choose h so that z=hλ falls inside the absolute stability region of your method. For a decaying system (λ<0), the stability limit is h < 2/|λ| for Euler and h < about 2.78/|λ| for RK4. Stability alone does not guarantee accuracy, however, so in practice an even smaller h is set from the required accuracy (local truncation error). Stability is the minimum condition for not blowing up; accuracy is how close the answer is — they are two different things.
A stiff equation is a system that mixes a very fast-decaying component (an eigenvalue with large |λ|) with slowly varying components. Even though the solution looks smooth, h must be kept below 2/|λ_max| just to solve the fast component stably, which makes the computation inefficient. Explicit methods such as Euler and RK4 have a finite stability region and are poor at stiff problems; the standard remedy is an A-stable implicit method like backward Euler, whose stability region covers the entire left half-plane. Set λ to -30 in this tool to feel the limit of explicit methods.
Real-World Applications
Structural vibration and transient response analysis: When the dynamic response of a building or machine is found by time integration, the higher-frequency modes have larger |λ| and impose a tight explicit stability condition h<2/|λ|. Refining the mesh raises the highest natural frequency and shrinks the required step size. In explicit impact and blast analysis (e.g. central difference), this stability limit appears as the "Courant condition", and the element size dictates the time step.
Chemical-reaction and combustion simulation: Solving many elementary reactions whose rate constants differ by orders of magnitude is a textbook stiff problem. Shrinking h to match the time constant of the fast reaction means an enormous number of steps to follow the slow reaction. In practice, instead of an explicit method like RK4, an implicit solver such as the backward differentiation formula (BDF) is the standard choice, freeing the computation from the stability constraint.
Numerical analysis of electrical circuits and control systems: In transient analysis of RC or op-amp circuits, small-time-constant parasitic capacitances produce large |λ|. SPICE-family simulators adopt the trapezoidal rule or backward Euler precisely to secure stability against this stiffness. When discretizing the state equations of a control system, whether the product z=hλ of the sampling period h and the eigenvalue falls inside the stability region is a design prerequisite.
Checking CAE solver settings: When an explicit dynamic analysis suddenly diverges, the time step has usually exceeded the stability limit. With an understanding of the stability region like this tool gives, you can diagnose that "it diverged after I refined the mesh" or "it diverged after I stiffened the material" simply means the eigenvalue λ grew and z=hλ left the region. Conversely, switching to an implicit method removes the need to shrink h for stability alone.
Common Misconceptions and Pitfalls
The biggest misconception is the belief that "a higher-order method is always more stable". RK4 is more accurate than Euler and its stability region does extend further along the real axis (-2.78 versus -2), but "higher order" does not mean "unconditionally stable". No matter how far you raise the order of an explicit Runge-Kutta method, the stability region stays finite and never covers the whole left half-plane. For a stiff equation, raising the order does not solve the essential problem — you must switch to an A-stable implicit method. Many "I used RK4 and it still diverges" questions stem from this misunderstanding.
Next, assuming that "stable means correct". The condition |R(z)| ≤ 1 only means "the numerical solution does not diverge"; it has nothing to do with accuracy. For instance, when z is near the edge of the stability region, |R(z)| is close to 1, and a component that should decay quickly lingers, or decays in an oscillatory way. It is not unusual for a stable solution to look nothing like the exact one. Even after passing a stability check, always halve the step size separately and confirm the solution does not change (that it has converged).
Finally, the misconception that "the test equation y'=λy is too simple to be relevant to real problems". It looks like a toy linear scalar equation, but when an actual coupled nonlinear ODE is linearized around an equilibrium point, each eigenvalue of the Jacobian matrix plays the role of λ. So stability analysis applies directly to real problems in the form "does z=hλ_i fall inside the stability region for every Jacobian eigenvalue λ_i?". The eigenvalue with the largest magnitude governs the step size — that is the essence of stiffness, and the reason the test-equation stability theory is genuinely useful in practice.
How to Use
Enter the number of eigenvalues (lamNum, range 1–5) and specify the real part (lamRange, typically –10 to 0 for stable systems)
Set the step size h in hRange (0.001 to 0.1) and number of steps hNum to define integration duration
Choose integration method from teRange: Forward Euler (order 1), RK2 (order 2), or RK4 (order 4)
Click Simulate to compute z = hλ, plot the stability region boundary, and overlay your operating point
Review |R(z)| amplification factor, stability verdict (stable/unstable), and final value error
Worked Example
Stiff chemical kinetics: λ = –50 s⁻¹, h = 0.01 s, 100 steps over 1 second, initial condition y₀ = 1. Forward Euler gives z = –0.5 (stable region), |R(z)| = 0.5, final value 0.0066, exact 0.0067, error 0.15%. RK4 with same h yields tighter damping. Increasing h to 0.05 places Euler at z = –2.5 (outside stability region), causing divergence to –32.8 versus exact –0.0067, error 489%.
Practical Notes
Forward Euler stable only for |z| ≤ 1; use h ≤ 0.02 with λ = –50 to stay safe in circuit transient analysis
RK4 extends stability to |z| ≤ 2.78; ideal for ODE systems in structural dynamics where λ ≈ –100 rad/s
Stiff systems (condition number >1000) require implicit methods; this simulator shows why explicit RK fails beyond critical h
Aerospace heat transfer: λ = –0.1 (slow cooling) allows large h ≈ 0.5, but λ = –200 (fast transient) demands h < 0.01 for accuracy