Math::RungeKutta.pm - Integrating Systems of Differential Equations
use Math::RungeKutta; # When working on data in an array ... sub dydt { my ($t, @y) = @_; # the derivative function my @dydt; ... ; return @dydt; } @y = @initial_y; $t=0; $dt=0.4; # the initial conditions # For automatic timestep adjustment ... while ($t < $tfinal) { ($t, $dt, @y) = &rk4_auto(\@y, \&dydt, $t, $dt, 0.00001); &display($t, @y); } # Or, for fixed timesteps ... while ($t < $tfinal) { ($t, @y) = &rk4(\@y, \&dydt, $t, $dt); # Merson's 4th-order method &display($t, @y); } # alternatively, though not so accurate ... ($t, @y) = &rk2(\@y, \&dydt, $t, $dt); # Heun's 2nd-order method # Or, working on data in a hash... sub dydt { my ($t, %y) = @_; # the derivative function my %dydt; ... ; return %dydt; } %y = %initial_y; $t=0; $dt=0.4; # the initial conditions # For automatic timestep adjustment on hashes ... while ($t < $tfinal) { ($t, $dt, %y) = &rk4_auto(\%y, \%dydt, $t, $dt, 0.00001); &display($t, %y); } # Or, for fixed timesteps on hashes ... while ($t < $tfinal) { ($t, %y) = &rk4(\%y, \%dydt, $t, $dt); # Merson's 4th-order method &display($t, %y); } # alternatively, though not so accurate ... ($t, %y) = &rk2(\%y, \%dydt, $t, $dt); # Heun's 2nd-order method # or, also available but not exported by default ... import qw(:ALL); ($t, @y) = &rk4_classical(\@y, \&dydt, $t, $dt); # Runge-Kutta 4th-order ($t, @y) = &rk4_ralston(\@y, \&dydt, $t, $dt); # Ralston's 4th-order # or similarly for data in hashes.
RungeKutta.pm offers algorithms for the numerical integration of simultaneous differential equations of the form
dY/dt = F(t,Y)
where Y is an array of variables whose initial values Y(0) are known, and F is a function known from the dynamics of the problem.
The Runge-Kutta methods all involve evaluating the derivative function F(t,Y) more than once, at various points within the timestep, and combining the results to reach an accurate answer for the Y(t+dt). This module only uses explicit Runge-Kutta methods; the implicit methods involve, at each timestep, solving a set of simultaneous equations involving both Y(t) and F(t,Y), and this is generally intractable.
Three main algorithms are offered. rk2 is Heun's 2nd-order Runge-Kutta algorithm, which is relatively imprecise, but does have a large range of stability which might be useful in some problems. rk4 is Merson's 4th-order Runge-Kutta algorithm, which should be the normal choice in situations where the step-size must be specified. rk4_auto uses the step-doubling method to adjust the step-size of rk4 automatically to achieve a specified precision; this saves much fiddling around trying to choose a good step-size, and can also save CPU time by automatically increasing the step-size when the solution is changing only slowly.
Perl is not the right language for high-end numerical integration like global weather simulation, colliding galaxies and so on (if you need something like this you could check out xmds). But as Gear says, "Many equations that are solved on digital computers can be classified as trivial by the fact that even with an inefficient method of solution, little computer time is used. Economics then dictates that the best method is the one that minimises the human time of preparation of the program."
This module has been designed to be robust and easy to use, and should be helpful in solving systems of differential equations which arise within a Perl context, such as economic, financial, demographic or ecological modelling, mechanical or process dynamics, etc.
Version 1.07
where the arguments are: \@y a reference to the array of initial values of variables, \%y a reference to the hash of initial values of variables, \&dydt a reference to the function calculating the derivatives, $t the initial time, $dt the timestep.
The algorithm used is that derived by Ralston, which uses Lotkin's bound on the derivatives, and minimises the solution error (gamma=3/4). It is also known as the Heun method, though unfortunately several other methods are also known under this name. Two function evaluations are needed per timestep, and the remaining error is in the 3rd and higher order terms.
rk2 returns ($t, @y) where $t and @y are now the new values at the completion of the timestep, or it returns ($t, %y) if called with the data in a hashref.
The arguments are the same as in rk2.
The algorithm used is that developed by Merson, which performs five function evaluations per timestep. For the same timestep, rk4 is much more accurate than rk4_classical, so the extra function evaluation is well worthwhile.
rk4 returns ($t, @y) where $t and @y are now the new values at the completion of the timestep.
In the I>epsilon> form the arguments are: \@y a reference to the array of initial values of variables or \%y a reference to the hash of initial values of variables, \&dydt a reference to the function calculating the derivatives, $t the initial time, $dt the initial timestep, $epsilon the errors per step will be about $epsilon*$ymax
In the errors form the last argument is: \@errors a reference to an array of maximum permissible errors, or \%errors a reference to a hash, accordingly.
The first $epsilon calling form is useful when all the elements of @y are in the same units and have the same typical size (e.g. y[10] is population aged 10-11 years, y[25] is population aged 25-26 years). The default value of the 4th argument is $epsilon = 0.00001.
The second errors form is useful otherwise (e.g. $y[1] is gross national product, $y[2] is interest rate, or $y{'gross national product'} and $y{'interest rate'} accordingly. In this calling form, the permissible errors are specified in absolute size for each variable; they won't get scaled at all.
rk4_auto adjusts the timestep automatically to give the required precision. It does this by trying one full-timestep, then two half-timesteps, and comparing the results. (Merson's method, as used by rk4, was devised to be able to give an estimate of the remaining local error; for the record, it is 0.2*($ynp1[i]-$eta4[i]) in each term. rk4_auto does not exploit this feature because it only works for linear dydt functions of the form Ay + bt.)
rk4_auto needs 14 function evaluations per double-timestep, and it has to re-do 13 of those every time it adjusts the timestep.
rk4_auto returns ($t, $dt, @y) where $t, $dt and @y are now the new values at the completion of the timestep, or ($t, $dt, %y) accordingly.
rk4_auto performs a double timestep within $dt, and returns the final values; the values as they were at the midpoint do not normally get returned. However, if you want to draw a nice smooth graph, or to update a nice smoothly-moving display, those values as they were at the midpoint would be useful to you. Therefore, rk4_auto_midpoint provides a way of retrieving them.
Note that you must call rk4_auto first, which returns the values at time $t+$dt, then rk4_auto_midpoint subsequently, which returns the values at $t+$dt/2, in other words you get the two sets of values out of their chronological order. Sorry about this. For example,
while ($t<$tfinal) { ($t, $dt, @y) = &rk4_auto(\@y, \&dydt, $t, $dt, $epsilon); ($t_midpoint, @y_midpoint) = &rk4_auto_midpoint(); &update_display($t_midpoint, @y_midpoint); &update_display($t, @y); }
rk4_auto_midpoint returns ($t, @y) where $t and @y were the values at the midpoint of the previous call to rk4_auto; or ($t, %y) accordingly.
You will pass this subroutine by reference as the second argument to rk2, rk4 and rk4_auto. The name doesn't matter of course. It must expect the following arguments: $t the time (in case the equations are time-dependent), @y the array of values of variables or %y the hash of values of variables.
It must return an array (or hash, accordingly) of the derivatives of the variables with respect to time.
The following routines are not exported by default, but are exported under the ALL tag, so if you need them you should:
import Math::RungeKutta qw(:ALL);
The arguments and the return values are the same as in rk2 and rk4.
The algorithm used is the classic, elegant, 4th-order Runge-Kutta method, using four function evaluations per timestep: k0 = dt * F(y(n)) k1 = dt * F(y(n) + 0.5*k0) k2 = dt * F(y(n) + 0.5*k1) k3 = dt * F(y(n) + k2) y(n+1) = y(n) + (k0 + 2*k1 + 2*k2 + k3) / 6
The algorithm used is that developed by Ralston, which optimises rk4_classical to minimise the error bound on each timestep. This module does not use it as the default 4th-order method rk4, because Merson's algorithm generates greater accuracy, which allows the timestep to be increased, which more than compensates for the extra function evaluation.
There are a couple of example scripts in the examples/ subdirectory of the build directory. You can use their code to help you get your first application going.
This script uses Term::Clui (arrow keys and Return, or q to quit) to offer a selection of algorithms, timesteps and error criteria for the integration of a simple sine/cosine wave around one complete cycle. This was the script used as a testbed during development.
This script uses the vt100 or xterm 'moveto' and 'reverse' sequences to display a little simulation of three-body gravity. It uses rk4_auto because a shorter timestep is needed when two bodies are close to each other. It also uses rk4_auto_midpoint to smooth the display. By changing the initial conditions you can experience how sensitively the outcome depends on them.
Alas, things can go wrong in numerical integration.
One of the most fundamental is instability. If you choose a timestep $dt much larger than time-constants implied in your derivative function &dydt, then the numerical solution will oscillate wildy, and bear no relation to the real behaviour of the equations. If this happens, choose a shorter $dt.
Some of the most difficult problems involve so-called stiff derivative functions. These arise when &dydt introduces a wide range of time-constants, from very short to long. In order to avoid instability, you will have to set $dt to correspond to the shortest time-constant; but this makes it impossibly slow to follow the evolution of the system over longer times. You should try to separate out the long-term part of the problem, by expressing the short-term process as the finding of some equilibrium, and then assume that that equilibrium is present and solve the long-term problem on its own.
Similarly, numerical integration doesn't enjoy problems where time-constants change suddenly, such as balls bouncing off hard surfaces, etc. You can often tackle these by intervening directly in the @y array between each timestep. For example, if $y[17] is the height of the ball above the floor, and $y[20] is the vertical component of the velocity, do something like
if ($y[17]<0.0) { $y[17]*=-0.9; $y[20]*=-0.9; }
and thus, again, let the numerical integration solve just the smooth part of the problem.
In the js/ subdirectory of the install directory there is RungeKutta.js, which is an exact translation of this Perl code into JavaScript. The function names and arguments are unchanged. Brief Synopsis:
js/
<SCRIPT type="text/javascript" src="RungeKutta.js"> </SCRIPT> <SCRIPT type="text/javascript"> var dydt = function (t, y) { // the derivative function var dydt_array = new Array(y.length); ... ; return dydt_array; } var y = new Array(); // For automatic timestep adjustment ... y = initial_y(); var t=0; var dt=0.4; // the initial conditions // Arrays of return vaules: var tmp_end = new Array(3); var tmp_mid = new Array(2); while (t < tfinal) { tmp_end = rk4_auto(y, dydt, t, dt, 0.00001); tmp_mid = rk4_auto_midpoint(); t=tmp_mid[0]; y=tmp_mid[1]; display(t, y); // e.g. could use wz_jsgraphics.js or SVG t=tmp_end[0]; dt=tmp_end[1]; y=tmp_end[2]; display(t, y); } // Or, for fixed timesteps ... y = post_ww2_y(); var t=1945; var dt=1; // start in 1945 var tmp = new Array(2); // Array of return values while (t <= 2100) { tmp = rk4(y, dydt, t, dt); // Merson's 4th-order method t=tmp[0]; y=tmp[1]; display(t, y); } </SCRIPT>
RungeKutta.js uses several global variables which all begin with the letters _rk_ so you should avoid introducing variables beginning with these characters.
_rk_
In the lua/ subdirectory of the install directory there is RungeKutta.lua, which is an exact translation of this Perl code into Lua. The function names and arguments are unchanged. Brief Synopsis:
lua/
local RK = require 'RungeKutta' function dydt(t, y) -- the derivative function -- y is the table of the values, dydt the table of the derivatives local dydt = {}; ... ; return dydt end y = initial_y(); t=0; dt=0.4; -- the initial conditions -- For automatic timestep adjustment ... while t < tfinal do t, dt, y = RK.rk4_auto(y, dydt, t, dt, 0.00001) display(t, y) end -- Or, for fixed timesteps ... while t < tfinal do t, y = RK.rk4(y, dydt, t, dt) -- Merson's 4th-order method display(t, y) end -- alternatively, though not so accurate ... t, y = RK.rk2(y, dydt, t, dt) -- Heun's 2nd-order method -- or, also available ... t, y = RK.rk4_classical(y, dydt, t, dt) -- Runge-Kutta 4th-order t, y = RK.rk4_ralston(y, dydt, t, dt) -- Ralston's 4th-order
Peter J Billam, http://www.pjb.com.au/comp/contact.html
On the Accuracy of Runge-Kutta's Method, M. Lotkin, MTAC, vol 5, pp 128-132, 1951
An Operational Method for the study of Integration Processes, R. H. Merson, Proceedings of a Symposium on Data Processing, Weapons Research Establishment, Salisbury, South Australia, 1957
Numerical Solution of Ordinary and Partial Differential Equations, L. Fox, Pergamon, 1962
A First Course in Numerical Analysis, A. Ralston, McGraw-Hill, 1965
Numerical Initial Value Problems in Ordinary Differential Equations, C. William Gear, Prentice-Hall, 1971
See also the scripts examples/sine-cosine and examples/three-body, http://www.pjb.com.au/, http://www.pjb.com.au/comp/, Math::WalshTransform, Math::Evol, Term::Clui, Crypt::Tea_JS, http://www.xmds.org/
To install Math::RungeKutta, copy and paste the appropriate command in to your terminal.
cpanm
cpanm Math::RungeKutta
CPAN shell
perl -MCPAN -e shell install Math::RungeKutta
For more information on module installation, please visit the detailed CPAN module installation guide.