view release on metacpan or  search on metacpan

RungeKutta.pm  view on Meta::CPAN

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 I<@y> array between each timestep. For example, if I<$y[17]>
is the height of the ball above the floor, and I<$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.

=head1 JAVASCRIPT

In the C<js/> subdirectory of the install directory there is I<RungeKutta.js>,
which is an exact translation of this Perl code into JavaScript.
The function names and arguments are unchanged.
Brief Synopsis:

 <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>

I<RungeKutta.js> uses several global variables
which all begin with the letters C<_rk_> so you should
avoid introducing variables beginning with these characters.

=head1 LUA

In the C<lua/> subdirectory of the install directory there is
I<RungeKutta.lua>, which is an exact translation of this Perl code into Lua.
The function names and arguments are unchanged.
Brief Synopsis:

 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

=head1 AUTHOR

Peter J Billam, http://www.pjb.com.au/comp/contact.html

=head1 REFERENCES

I<On the Accuracy of Runge-Kutta's Method>,
M. Lotkin, MTAC, vol 5, pp 128-132, 1951

I<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

I<Numerical Solution of Ordinary and Partial Differential Equations>,
L. Fox, Pergamon, 1962

I<A First Course in Numerical Analysis>, A. Ralston, McGraw-Hill, 1965

I<Numerical Initial Value Problems in Ordinary Differential Equations>,
C. William Gear, Prentice-Hall, 1971

=head1 SEE ALSO

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/

=cut