Mike Bostock’s page uses the Odex.js[1] library for running the differential equation. The library includes a predator-prey demo[2] of its own. I assumed Odex is the “gold standard”. I tried some quicker algorithms to see how close I could get.
First up, Euler’s method is the simplest numerical integration. Here’s a chart of how Euler’s method (green) compares to Odex (gray), both near the beginning of the simulation and after 100 seconds, with a tick size of 1/40th second:
We can see that Euler starts to diverge from the correct value 6 seconds into the simulation, and by the time we’re 100 seconds in, it’s pretty far off, not only in phase but also in amplitude. If I make the tick size smaller, Euler’s method will improve, but it will run slower.
What happens here is that in each peak and valley of the curve, dumb integration overshoots. {sketch} This is because it uses the slope from the left side (tangent) instead of the combined slope from the left and right (secant). We can’t calculate the slope using the secant without knowing the value on the right. And we can only calculate the value on the right by using the slope. So we’re kind of stuck.
Numerical integration has lots of issues like this, and I haven’t a lot of experience so I feel like I’m likely to get bitten. Despite that, I want to draw enough charts that I want to get an approximation very quickly than get the right answer slowly.
The next thing to try is one of the predictor-corrector[3] methods, which use Euler’s method as an approximation to construct the secant, and then use the secant for the slope. Here’s a chart of how Heun’s method[4] (red) compares to Odex (gray), near the beginning of the simulation, after 100 seconds, after 1000 seconds, and even after 10000 seconds:
Heun’s method is much more accurate than Euler’s method! It’s a little slower than Euler’s but both are much faster than Odex. I decided to run it for longer, up to 1000 and 10000 seconds:
Heun’s is slightly different from Odex around 1000 seconds in, but I think it’s close enough that I’m not going to worry about it. By 10000 seconds it’s quite different from Odex. I’m glad the amplitude doesn’t diverge like Euler’s does. And 1/40th of a second seems to work well.
The next thing I wanted to look at was the phase plots (predator vs prey instead of vs time). The solutions should be periodic. If I plot them over time, they should be exactly in the same place every time through the oscillation. Euler’s method is clearly not after 100 seconds:
Heun’s method is better, black line is at the beginning and red line after 10,000 seconds:
With Odex though, the red line after 10,000 seconds differs from the black line at the beginning:
For this particular problem, Heun’s method appears to be better than Odex. I wanted to investigate more, so I plotted the maximum value over time. Since it’s periodic, it should never go up. I already knew Euler’s method was diverging quickly. Within 100 seconds you can see every time it orbits, the amplitude has increased:
Here’s Heun’s method (red line) and Odex (black line), showing that both have increasing amplitude, but Heun’s method is more stable:
I think that’s enough for now. I know what I wanted to know: Heun’s is not only good enough, it’s probably a little better than Odex for this particular problem.