Stiff Systems: Adaptive Steps and Implicit Methods
Claude prompt
This page was generated by Claude Code from the following prompt (with additional prompting):
let's build a van-der-pol stiff system javascript demo with BS23 and implicit Euler, let's also add fixed-step RK4 and explicit euler to see things that blow up
Test problem: the van der Pol oscillator, \[ \ddot{x} = \mu\,(1 - x^2)\,\dot{x} - x, \qquad x(0) = 2,\quad \dot{x}(0) = 0. \] For small \(\mu\) it is a gentle oscillator. For large \(\mu\) it becomes stiff: the solution creeps along slowly for most of each cycle, then jumps almost instantaneously. The reference curve is BS23 with \(E_s = 10^{-9}\).
Fixed \(h\). Stable only while \(|1 + h\lambda| < 1\) for every eigenvalue \(\lambda\) of the Jacobian, i.e. \(h < 2/|\lambda|\) for real negative \(\lambda\).
Fixed \(h\). Fourth order, but still explicit: the stability limit is \(h < 2.785/|\lambda|\), barely better than Euler's.
Bogacki–Shampine 3(2), the method inside MATLAB's ode23, with the error estimate from an embedded second-order step. Still explicit: when a step crosses the stability limit the error estimate explodes and the step is rejected, so the method never blows up. It just crawls.
Each step solves a 2×2 nonlinear system for \(\vec{x}_{i+1}\) with Newton's method. Step halving gives the error estimate. First order, so it needs small steps for accuracy, but it is stable for any \(h\) when \(\mathrm{Re}\,\lambda < 0\), so the step size is set by accuracy alone.
A Rosenbrock method, what MATLAB's ode23s uses. Implicit, second order, with a built-in error estimate.
Trajectories
Black: reference solution. Colored: the selected method, with a dot at each step when there are few enough to see. The red bar at \(t_f\) is the error in each component.
Efficiency of the adaptive methods at the current \(\mu\): error at \(t_f\) against the number of evaluations of \(\vec{f}\), one point per value of \(E_s\). Down and to the left is better. The red circle is the current setting.