Integration of ODEs with Adaptive Step Size using Perl

I have a strong preference for Monte Carlo simulations over deterministic simulations, because the former allow for a much more simple and direct representation of the physics and are far less susceptible to simple mistakes.

I once worked on an experiment that had been taking data for over a year in a configuration that had a background which perfectly imitated the signal. The original design calculations contained a sign error that caused the cancellation of two terms that in fact added together. Monte Carlo simulation showed the large background, and this was confirmed by independent deterministic integration. Once he knew there was an error, the original designer was able to find the droppe sign in his calculation.

I know for a fact he was a careful physicist; it wasn't sloppiness that caused the problem, but the inherent difficulty of keeping the abstract representation correct.

While tools like Mathematica (one of my favorite applications) now save us from much error-prone algebraic tedium, deterministic simulations are still subject to representational errors, particularly due to simple slips like writing the force on something instead of the force of something. These sorts of mistakes manifest themselves as sign errors even when all the rest of the algebra is correct.

Despite all this, Monte Carlo simulations are not enough: they suffer from the inverse problem of deterministic simulations. The latter can be too abstract to admit of simple interpretation. The former can be too concrete, making it impossible to see the forest for the trees, or worse yet, the leaves.

Some of my current work involves an investigation of the relationship between classical thermodynamics and kinetic theory. This is still, a hundred years after Boltzmann's death, a thorny conceptual issue: just how far can the laws of mechanics take us toward classical thermodynamics? What is the best way to answer Loschmidt's objections to Boltzmann's justification of thermodynamic laws in terms of ergodic mechanics? There is a gap in our understanding that we paper over too easily, a gap that seems tantalizingly bridgeable.

Monte Carlo modeling allows us to investigate the interactions of individual particles in great detail, which can shed some light on these questions. But to really illuminate the issues we need to perform a parallel analysis in terms of equations of motion, and to solve those differential equations it is convenient to have a simple ODE solver lying around.

I've written ODE solvers in the past, but for various reasons none of them are still available to me, so I threw together something quick and dirty for the system I was working on in Perl. Perl isn't a natural language for numerical analysis, but I find when taking multiple approaches to a given problem it helps keep them distinct and compartmentalized if I use a different language for each approach. I was already using C++ for the Monte Carlo, and I'm really off FORTRAN right now despite it being my language of choice for more than a decade, so I decided Perl was the thing. It seemed like a good idea at the time.

I found John Leto's Math::ODE module (v0.03) on CPAN, and it worked nicely for the simple case I was interested in, which was nothing more than a linear harmonic oscillator. That case was a linearized version of the real system of interest (it always helps to fall back on something you can solve with pen and paper--computers produce numbers, not insight!) The real system is nonlinear and the equations of motion have singularities at the extrema.

Math::ODE v0.03 is a simple 2nd order Runge-Kutta solver, and for any finite step size under some circumstances the motion would come too close to a singularity and the solution would go off the rails. Clearly the solution was to add an adaptive step-size to the module, which is a pretty good thing to have for even more mundane applications of Runge-Kutta.

The Math::ODE v0.04 module is available here for download. It contains the step-size adapter as well as a few other mostly cosmetic modifications. I'm in the midst of communicating with the original author to see about getting it submitted to CPAN for the rest of the world to enjoy.