#!/usr/bin/env ruby
# Ruby/GSL implementation of GSL "odeiv-initval/test.c"
require("rbgsl")
require("./gsl_test2.rb")
include GSL::Test
include Math
rhs_linear = Proc.new { |t, y, f|
f[0] = 0.0
f[1] = y[0]
GSL::SUCCESS
}
jac_linear = Proc.new { |t, y, dfdy, dfdt|
dfdy.set(0, 0, 0.0)
dfdy.set(0, 1, 0.0)
dfdy.set(1, 0, 1.0)
dfdy.set(1, 1, 0.0)
dfdt[0] = 0.0
dfdt[1] = 0.0
GSL::SUCCESS
}
Rhs_func_lin = GSL::Odeiv::System.alloc(rhs_linear, jac_linear, 2)
rhs_sin = Proc.new { |t, y, f|
f[0] = -y[1]
f[1] = y[0]
GSL::SUCCESS
}
jac_sin = Proc.new { |t, y, dfdy, dfdt|
dfdy.set(0, 0, 0.0)
dfdy.set(0, 1, -1.0)
dfdy.set(1, 0, 1.0)
dfdy.set(1, 1, 0.0)
dfdt[0] = 0.0
dfdt[1] = 0.0
GSL::SUCCESS
}
Rhs_func_sin = GSL::Odeiv::System.alloc(rhs_sin, jac_sin, 2)
rhs_exp = Proc.new { |t, y, f|
f[0] = y[1]
f[1] = y[0]
GSL::SUCCESS
}
jac_exp = Proc.new { |t, y, dfdy, dfdt|
dfdy.set(0, 0, 0.0)
dfdy.set(0, 1, 1.0)
dfdy.set(1, 0, 1.0)
dfdy.set(1, 1, 0.0)
dfdt[0] = 0.0
dfdt[1] = 0.0
GSL::SUCCESS
}
Rhs_func_exp = GSL::Odeiv::System.alloc(rhs_exp, jac_exp, 2)
rhs_stiff = Proc.new { |t, y, f|
f[0] = 998.0 * y[0] + 1998.0 * y[1]
f[1] = -999.0 * y[0] - 1999.0 * y[1]
GSL::SUCCESS
}
jac_stiff = Proc.new { |t, y, dfdy, dfdt|
dfdy.set(0, 0, 998.0);
dfdy.set(0, 1, 1998.0);
dfdy.set(1, 0, -999.0);
dfdy.set(1, 1, -1999.0);
dfdt[0] = 0.0;
dfdt[1] = 0.0;
GSL::SUCCESS;
}
Rhs_func_stiff = GSL::Odeiv::System.alloc(rhs_stiff, jac_stiff, 2)
def test_stepper_linear(type, hstart, base_prec)
h = hstart*1.0
delmax = 0.0
count = 0
stepper = GSL::Odeiv::Step.alloc(type, 2)
y = GSL::Vector.alloc(1.0, 0.0)
yerr = GSL::Vector.alloc(2)
t = 0.0
s = 0
while t < 4.0
status = stepper.apply(t, h, y, yerr, Rhs_func_lin)
del = ((y[1] - (t + h))/y[1]).abs
delmax = MAX_DBL(del, delmax)
if del > (count + 1.0)*base_prec
printf(" LINEAR(%20.17g) %20.17g %20.17g %8.4g\n", t + h, y[1], t + h, del);
s += 1
end
count += 1
t += h
end
GSL::Test::test(s, "#{stepper.name}, linear [0,4], max relative error = #{delmax}")
end
def test_stepper_sin(type, hstart, base_prec)
h = hstart*1.0
delmax = 0.0
count = 0
stepper = GSL::Odeiv::Step.alloc(type, 2)
y = GSL::Vector.alloc(1.0, 0.0)
yerr = GSL::Vector.alloc(2)
t = 0.0
s = 0
while t < M_PI
sin_th = sin(t + h)
status = stepper.apply(t, h, y, yerr, Rhs_func_sin)
del = ((y[1] - sin_th)/sin_th).abs
delmax = GSL::MAX_DBL(del, delmax)
if t < 0.5*M_PI
stat = (del > (count+1.0)*base_prec) ? 1 : 0
elsif t < 0.7*M_PI
stat = (del > 1.0e+04*base_prec) ? 1 : 0
elsif t < 0.9*M_PI
stat = (del > 1.0e+06*base_prec) ? 1 : 0
else
stat = (del > 1.0e+09*base_prec) ? 1 : 0
end
if stat != 0
printf(" SIN(%22.18g) %22.18g %22.18g %10.6g\n", t + h, y[1], sin_th, del);
end
s += stat
count += 1
t += h
end
if delmax > 1.0e+09*base_prec
s += 1
printf(" SIN(0 .. M_PI) delmax = %g\n", delmax)
end
GSL::Test::test(s, "#{stepper.name}, sine [0,pi], max relative error = #{delmax}")
delmax = 0.0
while t < 3*M_PI
status = stepper.apply(t, h, y, yerr, Rhs_func_sin)
del = (y[1] - sin(t)).abs
delmax = GSL::MAX_DBL(del, delmax)
count += 1
t += h
end
if del > count * 2.0 * base_prec
s += 1
printf(" SIN(%22.18g) %22.18g %22.18g %10.6g\n", t + h, y[1], sin(t), del)
end
test(s, "#{stepper.name}, sin [pi,3*pi], max absolute error = #{delmax}")
end
def test_stepper_exp(type, hstart, base_prec)
s = 0
y = GSL::Vector.alloc(1.0, 1.0)
yerr = GSL::Vector.alloc(2)
delmax = 0.0
count = 0
h = hstart*1.0
stepper = GSL::Odeiv::Step.alloc(type, 2)
t = 0.0
while t < 5.0
ex = exp(t + h)
status = stepper.apply(t, h, y, yerr, Rhs_func_exp)
del = ((y[1] - ex)/y[1]).abs
delmax = GSL::MAX_DBL(del, delmax)
if del > (count+1.0)*2.0*base_prec
printf(" EXP(%20.17g) %20.17g %20.17g %8.4g\n", t + h, y[1], ex, del);
s += 1
end
count += 1
t += h
end
GSL::Test::test(s, "#{stepper.name}, exponential [0,5], max relative error = #{delmax}")
end
def test_stepper_stiff(type, hstart, base_prec)
h = hstart*1.0
s = 0
delmax = 0.0
count = 0
stepper = GSL::Odeiv::Step.alloc(type, 2)
y = GSL::Vector.alloc(1.0, 0.0)
yerr = GSL::Vector.alloc(2)
t = 0.0
while t < 5.0
status = stepper.apply(t, h, y, yerr, Rhs_func_stiff)
if t > 0.04
arg = t + h
e1 = exp(-arg)
e2 = exp(-1000.0*arg)
u = 2.0*e1 - e2
del = ((y[0] - u)/y[0]).abs
delmax = GSL::MAX_DBL(del, delmax)
if del > (count + 1.0)*100.0*base_prec
printf(" STIFF(%20.17g) %20.17g %20.17g %8.4g\n", arg, y[0], u, del)
s += 1
end
end
count += 1
t += h
end
GSL::Test::test(s, "#{stepper.name}, stiff [0,5], max relative error = #{delmax}")
end
def test_stepper_err(type, hstart, base_prec)
h = hstart*1.0
y = GSL::Vector.alloc(1.0, 0.0)
yerr = GSL::Vector.alloc(2)
delmax = 0.0
errmax = 0.0
count = 0
s = 0
stepper = GSL::Odeiv::Step.alloc(type, 2)
t = 0.0
while t < M_PI
y1_t = y[1]
dy_exp = cos(t)*sin(h)-2*sin(t)*pow(sin(h/2),2.0);
status = stepper.apply(t, h, y, yerr, Rhs_func_sin)
dy_t = y[1] - y1_t;
del = (dy_t - dy_exp).abs
if t > 0.1 and t < 0.2
stat = (del > 10.0*yerr[1].abs + GSL::DBL_EPSILON*y[1].abs) ? 1 : 0
if stat != 0
delmax = del
errmax = yerr[1].abs
printf("SIN(%.18e) %.5e %.5e %e %e %e\n", t + h, y[1], dy_t, dy_exp, del, yerr[1]);
s += stat;
break;
end
end
count += 1
t += h
end
test(s, "#{stepper.name}, sine [0,pi], accurary of estimate error = #{delmax} vs #{errmax}")
end
def test_evolve_system_flat(step, sys, t0, t1, hstart, y, yfin, err_target, desc)
s = 0
h = hstart*1.0
t = t0*1.0
e = GSL::Odeiv::Evolve.alloc(sys.dimension)
while t < t1
t, h, status = e.apply(nil, step, sys, t, t1, h, y)
end
frac = ((y[1] - yfin[1])/yfin[1]).abs + ((y[0] - yfin[0])/yfin[0]).abs
if frac > 2.0*e.count*err_target
printf("FLAT t = %.5e y0 = %g y1= %g\n", t, y[0], y[1]);
s += 1
end
GSL::Test::test(s, "#{step.name}, #{desc}, evolve, no control, max relative error = #{frac}")
end
def test_evolve_system(type, sys, t0, t1, hstart, y, yfin, err_target, desc)
s = 0
t = t0*1.0
h = hstart*1.0
step = GSL::Odeiv::Step.alloc(type, sys.dimension)
c = GSL::Odeiv::Control.standard_alloc(0.0, err_target, 1.0, 1.0)
e = GSL::Odeiv::Evolve.alloc(sys.dimension)
while t < t1
t, h, status = e.apply(c, step, sys, t, t1, h, y)
end
frac = ((y[1] - yfin[1])/yfin[1]).abs + ((y[0] - yfin[0])/yfin[0]).abs
if frac > 2.0*e.count*err_target
printf("SYS t = %.5e h = %g y0 = %g y1= %g\n", t, h, y[0], y[1]);
s += 1
end
GSL::Test::test(s, "#{step.name}, #{desc}, evolve, standard control, relative error = #{frac}")
end
def test_evolve_exp(type, hstart, err)
h = hstart*1.0
y = GSL::Vector.alloc(1.0, 1.0)
eee = exp(5.0)
yfin = GSL::Vector.alloc(eee, eee)
test_evolve_system(type, Rhs_func_exp, 0.0, 5.0, h, y, yfin, err, "exp [0,5]")
end
def test_evolve_sin(type, hstart, err)
h = hstart*1.0
y = GSL::Vector.alloc(1.0, 0.0)
yfin = GSL::Vector.alloc(cos(2.0), sin(2.0))
test_evolve_system(type, Rhs_func_sin, 0.0, 2.0, h, y, yfin, err, "sin [0,2]")
end
def test_evolve_stiff1(type, hstart, err)
h = hstart*1.0
y = GSL::Vector.alloc(1.0, 0.0)
arg = 1.0
e1 = exp(-arg)
e2 = exp(-1000.0*arg)
yfin = GSL::Vector.alloc(2.0*e1 - e2, -e1 + e2)
test_evolve_system(type, Rhs_func_stiff, 0.0, 1.0, h, y, yfin, err, "stiff [0,1]")
end
def test_evolve_stiff5(type, hstart, err)
h = hstart*1.0
y = GSL::Vector.alloc(1.0, 0.0)
arg = 5.0
e1 = exp(-arg)
e2 = exp(-1000.0*arg)
yfin = GSL::Vector.alloc(2.0*e1 - e2, -e1 + e2)
less test_evolve_system(type, Rhs_func_stiff, 0.0, 5.0, h, y, yfin, err, "stiff [0,5]")
end
ptypes = [
{"type" => "rk2",
"h" => 1e-3},
{"type" => "rk2imp",
"h" => 1e-3},
{"type" => "rk4",
"h" => 1e-3},
{"type" => "rk4imp",
"h" => 1e-3},
{"type" => "rkf45",
"h" => 1e-3},
{"type" => "rk8pd",
"h" => 1e-3},
{"type" => "rkck",
"h" => 1e-3},
{"type" => "bsimp",
"h" => 1e-3},
{"type" => "gear1",
"h" => 1e-3},
{"type" => "gear2",
"h" => 1e-3}
]
GSL::IEEE::env_setup()
ptypes.each do |hash|
test_stepper_err(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
end
ptypes.each do |hash|
# test_stepper_linear(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
# test_stepper_sin(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
test_stepper_sin(hash["type"], hash["h"]/10, 1e-8)
# test_stepper_exp(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
# test_stepper_stiff(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
end
ptypes.each do |hash|
# test_evolve_exp(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
test_evolve_sin(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
# test_evolve_stiff1(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
# test_evolve_stiff5(hash["type"], hash["h"], GSL::SQRT_DBL_EPSILON)
end
syntax highlighted by Code2HTML, v. 0.9.1