view release on metacpan or  search on metacpan

lua/RungeKutta.lua  view on Meta::CPAN

-- RungeKutta.lua
-- ----------------------------------------------------------------- --
--      This Lua5 module is Copyright (c) 2010, Peter J Billam       --
--                        www.pjb.com.au                             --
--                                                                   --
--   This module is free software; you can redistribute it and/or    --
--          modify it under the same terms as Lua5 itself.           --
-- ----------------------------------------------------------------- --
local M = {} -- public interface
M.Version = '1.07'
M.VersionDate = '20aug2010'

-- Example usage:
-- local RK = require 'RungeKutta'
-- RK.rk4()

--------------------- infrastructure ----------------------
local function arr2txt(...) -- neat printing of arrays for debug use
	local txt = {}
	for e in ... do txt[#txt+1] = string.format('%g',e) end
	return table.concat(txt,' ') .. "\n"
end
local function warn(str)
	io.stderr:write(str,'\n')
end
local function die(str)
	io.stderr:write(str,'\n')
	os.exit(1)
end
local flag = false
local a
local b
local function gaussn(standdev)
	-- returns normal distribution around 0.0 by the Box-Muller rules
	if not flag then
		a = math.sqrt(-2.0 * math.log(math.random()))
		b = 6.28318531 * math.random()
		flag = true
		return standdev * a * math.sin(b)
	else
		flag = false
		return standdev * a * math.cos(b)
	end
end
-------------------------------------------------------

function M.rk2(yn, dydt, t, dt)
	if type(yn) ~= 'table' then
		warn("RungeKutta.rk2: 1st arg must be an table\n")
		return false
	end
	if type(dydt) ~= 'function' then
		warn("RungeKutta.rk2: 2nd arg must be a function\n")
		return false
	end

	local gamma = .75;  -- Ralston's minimisation of error bounds
	local alpha = 0.5/gamma; local beta = 1.0-gamma;
	local alphadt=alpha*dt; local betadt=beta*dt; local gammadt=gamma*dt;
	local ny = #yn;
	local ynp1 = {}
	local dydtn = {}
	local ynpalpha = {}  -- Gear calls this q
	local dydtnpalpha = {}
	dydtn = dydt(t, yn);
	-- for i=1, ny do
	for i in pairs(yn) do
		ynpalpha[i] = yn[i] + alphadt*dydtn[i];
	end
	dydtnpalpha = dydt(t+alphadt, ynpalpha);
	for i in pairs(yn) do
		ynp1[i] = yn[i]+betadt*dydtn[i]+gammadt*dydtnpalpha[i];
	end
	return t+dt, ynp1

lua/RungeKutta.lua  view on Meta::CPAN

	for i in pairs(yn) do k3[i] = dt * k3[i] end

	local ynp1 = {}
	for i in pairs(yn) do
		ynp1[i] = yn[i] + 0.17476028*k0[i]
		 - 0.55148053*k1[i] + 1.20553547*k2[i] + 0.17118478*k3[i]
	end
	return t+dt, ynp1
end
function M.rk4_classical(yn, dydt, t, dt)
	if (type(yn) ~= 'table') then
		warn("RungeKutta.rk4_classical: 1st arg must be arrayref\n")
		return false
	end
	if (type(dydt) ~= 'function') then
		warn("RungeKutta.rk4_classical: 2nd arg must be subroutine ref\n")
		return false
	end
	local ny = #yn; local i;

	-- The Classical 4th-order Runge-Kutta Method, see Gear p35

	local k0 = dydt(t, yn)
	for i in pairs(yn) do k0[i] = dt * k0[i] end

	local eta1 = {}
	for i in pairs(yn) do eta1[i] = yn[i] + 0.5*k0[i] end
	local k1 = dydt(t+0.5*dt, eta1)
	for i in pairs(yn) do k1[i] = dt * k1[i] end

	local eta2 = {}
	for i in pairs(yn) do eta2[i] = yn[i] + 0.5*k1[i] end
	local k2 = dydt(t+0.5*dt, eta2)
	for i in pairs(yn) do k2[i] = dt * k2[i] end

	local eta3 = {}
	for i in pairs(yn) do eta3[i] = yn[i] + k2[i] end
	local k3 = dydt(t+dt, eta3)
	for i in pairs(yn) do k3[i] = dt * k3[i] end

	local ynp1 = {}
	for i in pairs(yn) do
		ynp1[i] = yn[i] + (k0[i] + 2.0*k1[i] + 2.0*k2[i] + k3[i]) / 6.0;
	end
	return t+dt, ynp1
end


return M

--[[

=pod

=head1 NAME

RungeKutta.lua - Integrating Systems of Differential Equations

=head1 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
   -- the table can be an array (1...n), or a dictionary; whichever,
   -- the same indices must be used for the return table: dydt
   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 DESCRIPTION

RungeKutta.lua 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.  I<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.  I<rk4>
is Merson's 4th-order Runge-Kutta algorithm, which should be the normal
choice in situations where the step-size must be specified.  I<rk4_auto>
uses the step-doubling method to adjust the step-size of I<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.

This module is the translation into I<Lua> of the I<Perl> CPAN module
Math::RungeKutta, and comes in its C<./lua> subdirectory.
There also exists a translation into I<JavaScript>
which comes in its C<./js> subdirectory.
The calling-interfaces are identical in all three versions.

This module has been designed to be robust and easy to use, and should
be helpful in solving systems of differential equations which arise