#!/usr/bin/lua
-- ----------------------------------------------------------------- --
-- This Lua5 script is Copyright (c) 2010, Peter J Billam --
-- www.pjb.com.au --
-- --
-- This script is free software; you can redistribute it and/or --
-- modify it under the same terms as Lua5 itself. --
-- ----------------------------------------------------------------- --
local Version = '1.0 for Lua5'
local VersionDate = '31jul2010';
local Synopsis = [[
program_name [options] [filenames]
]]
local RK = require 'RungeKutta'
--------------------------- infrastructure -----------------------
local eps = .000000001
function equal(x, y) -- unused here
if #x ~= #y then return false end
local i; for i in pairs(x) do
if math.abs(x[i]-y[i]) > eps then return false end
end
return true
end
-- use Test::Simple tests => 6;
local Test = 12 ; local i_test = 0; local Failed = 0;
function ok(b,s)
i_test = i_test + 1
if b then
io.write('ok '..i_test..' - '..s.."\n")
else
io.write('not ok '..i_test..' - '..s.."\n")
Failed = Failed + 1
end
end
------------------------------------------------------------------
local func_evals = 0
local i_test = 0
local n_failed = 0
local n_passed = 0
function dydt(t, y)
local dydt_f = {}
dydt_f[1] = y[2]
dydt_f[2] = 0.0 - y[1]
func_evals = func_evals + 1
return dydt_f
end
function dydt_hash(t, ya)
local i
local dydt_f = {}
dydt_f["x"] = 0.0 + ya["y"]
dydt_f["y"] = 0.0 - ya["x"]
func_evals = func_evals + 1
return dydt_f
end
twopi = 2.0 * 3.141592653589
algorithms = {RK.rk2, RK.rk4, RK.rk4_classical, RK.rk4_ralston, epsilon, errors}
algnames = {'rk2', 'rk4', 'rk4_classical', 'rk4_ralston'}
passmark0 = { 0.2, .0004, .0015, .0015, .0001, .0003 }
passmark1 = { 0.04, .00004, .0006, .0006, .00001, .00001 }
for i=1,4 do
algorithm = algorithms[i]
i_test = i_test + 1
n = 16
dt= twopi / n
y = {0,1}; t=0
for j=1, n do t, y = algorithm(y, dydt, t, dt) end
local err0 = math.abs(y[1]);
local err1 = math.abs(y[2]-1.0)
ok(err0 < passmark0[i] and err1 < passmark1[i], algnames[i]..' with array')
y = {x=0, y=1}; t=0
for j=1, n do t, y = algorithm(y, dydt_hash, t, dt) end
local err0 = math.abs(y['x']);
local err1 = math.abs(y['y']-1.0)
ok(err0 < passmark0[i] and err1 < passmark1[i],algnames[i]..' with hash')
end
algorithm = RK.rk4_auto
local t_midpoint; local y_midpoint
for k, mode in ipairs({'epsilon','errors'}) do
local skip_to_next = false
i_test = i_test + 1
local i = 0
local epsilon
--- array ---
if mode == 'epsilon' then epsilon = .0001
else errors = {.01, .0001}; epsilon = errors;
end
y = {0,1}; t = 0; dt = 0.1
func_evals = 0
while t+dt < twopi do
i = i + 1
t, dt, y = RK.rk4_auto(y, dydt, t, dt, epsilon )
t_midpoint, y_midpoint = RK.rk4_auto_midpoint()
if (func_evals > 500) then ok(false,
'rk4_auto with array and '..mode..', '..func_evals..' func evals')
skip_to_next = true ; break
end
end
if not skip_to_next then
i = i + 1; dt = twopi-t;
t, y = RK.rk4(y, dydt, t, dt)
local err0 = math.abs(y[1]);
local err1 = math.abs(y[2]-1.0)
ok(err0 < passmark0[i_test] and err1 < passmark1[i_test],
'rk4_auto with array and '..mode)
end
--- hash ---
if mode == 'epsilon' then epsilon = .0001
else errors = {x=.01, y=.0001}; epsilon = errors;
end
y = {x=0,y=1}; t = 0; dt = 0.1
func_evals = 0
while t+dt < twopi do
i = i + 1
t, dt, y = RK.rk4_auto(y, dydt_hash, t, dt, epsilon )
t_midpoint, y_midpoint = RK.rk4_auto_midpoint()
if (func_evals > 500) then
ok(false,
'rk4_auto with hash and '..mode..', '..func_evals..' func evals')
skip_to_next = true ; break
end
end
if not skip_to_next then
i = i + 1; dt = twopi-t;
t, y = RK.rk4(y, dydt_hash, t, dt)
local err0 = math.abs(y["x"]);
local err1 = math.abs(y["y"]-1.0)
ok(err0 < passmark0[i_test] and err1 < passmark1[i_test],
'rk4_auto with hash and '..mode)
end
end
--[[
__END__
=pod
=head1 NAME
test.lua - Lua script to test RungeKutta.lua
=head1 SYNOPSIS
perl test.pl
=head1 DESCRIPTION
This script tests Math::RungeKutta.pm
=head1 AUTHOR
Peter J Billam http://www.pjb.com.au/comp/contact.html
=head1 SEE ALSO
Math::RungeKutta.pm , http://www.pjb.com.au/ , perl(1).
=cut
]]