| -rw-r--r-- | benchmarks/integration/qag-bench.c | 47 | ||||
| -rw-r--r-- | benchmarks/integration/qag-bench.lua | 28 | ||||
| -rw-r--r-- | benchmarks/odeiv/odeiv-benchmark.c (renamed from ode-evolve-benchmark.c) | 0 | ||||
| -rw-r--r-- | benchmarks/odeiv/odeiv-benchmark.lua (renamed from ode-evolve-benchmark.lua) | 2 | ||||
| -rw-r--r-- | benchmarks/odeiv/odeiv-vec-benchmark.lua (renamed from ode-evolve-benchmark-vec.lua) | 2 | ||||
| -rw-r--r-- | examples/ode-jit.lua | 129 | ||||
| -rw-r--r-- | examples/qag-bessel.lua | 17 | ||||
| -rw-r--r-- | examples/qng-bessel.lua | 27 | ||||
| -rw-r--r-- | num/gauss-kronrod-x-wgs.lua.in | 1 | ||||
| -rw-r--r-- | num/ode-defs.lua.in (renamed from ode-defs.lua.in) | 0 | ||||
| -rw-r--r-- | num/qag.lua.in | 21 | ||||
| -rw-r--r-- | num/qng.lua | 4 | ||||
| -rw-r--r-- | num/qng.lua.in | 4 | ||||
| -rw-r--r-- | num/rk4.lua.in (renamed from rk4.lua.in) | 2 | ||||
| -rw-r--r-- | num/rkf45.lua.in (renamed from rkf45.lua.in) | 33 | ||||
| -rw-r--r-- | num/rkf45vec.lua.in (renamed from rkf45vec.lua.in) | 0 | ||||
| -rw-r--r-- | ode-test-vec.lua | 104 | ||||
| -rw-r--r-- | ode-test.lua | 97 | ||||
| -rw-r--r-- | rk4-ref-evolve.c | 70 | ||||
| -rw-r--r-- | rk4-ref.c | 60 | ||||
| -rw-r--r-- | rk4.lua | 118 | ||||
| -rw-r--r-- | rkf45vec-idamax.lua.in | 181 |
diff --git a/benchmarks/integration/qag-bench.c b/benchmarks/integration/qag-bench.c new file mode 100644 index 00000000..4bfa3939 --- /dev/null +++ b/benchmarks/integration/qag-bench.c @@ -0,0 +1,47 @@ +#include <stdio.h> +#include <math.h> +#include <gsl/gsl_integration.h> + +struct bessel_param { + double x; + int n; +}; + +double f (double t, void * params) { + struct bessel_param *p = (struct bessel_param *) params; + double x = p->x; + int n = p->n; + return cos(n * t - x * sin(t)); +} + +int +main (void) +{ + gsl_integration_workspace * ws = gsl_integration_workspace_alloc (1000); + + struct bessel_param param = {0.0, 12}; + double result, error; + double xold = -100, xsmp = 10; + int k; + + gsl_function F; + F.function = &f; + F.params = ¶m; + + for (k = 0; k < 4096 * 8; k++) + { + param.x = (k * 30 * M_PI) / (4096 * 8); + + gsl_integration_qag (&F, 0, M_PI, 1e-6, 1e-4, 1000, GSL_INTEG_GAUSS21, ws, &result, &error); + + if (param.x - xold > xsmp) + { + xold = param.x; + printf ("%.18f %.18f\n", param.x, result); + } + } + + gsl_integration_workspace_free (ws); + + return 0; +} diff --git a/benchmarks/integration/qag-bench.lua b/benchmarks/integration/qag-bench.lua new file mode 100644 index 00000000..0cb7da3f --- /dev/null +++ b/benchmarks/integration/qag-bench.lua @@ -0,0 +1,28 @@ +local template = require 'template' +local q = template.load('num/qag.lua.in', {limit=1000, order=21}) + +local sin, cos, pi = math.sin, math.cos, math.pi +local epsabs, epsrel = 1e-6, 0.01 + +function bessel_gen(n) + local xs + local fint = function(t) return cos(n*t - xs*sin(t)) end + return function(x) + xs = x + return q(fint, 0, pi, epsabs, epsrel) + end +end + +local J12 = bessel_gen(12) + +-- p = graph.fxplot(J12, 0, 30*pi) + +local xold, xsmp = -100, 10 +for k = 1, 4096*8 do + local x = (k-1) * 30 * math.pi / (4096*8) + local y = J12(x); + if x - xold > xsmp then + print(x, y) + xold = x + end +end diff --git a/ode-evolve-benchmark.c b/benchmarks/odeiv/odeiv-benchmark.c index 2b9c92ed..2b9c92ed 100644 --- a/ode-evolve-benchmark.c +++ b/benchmarks/odeiv/odeiv-benchmark.c diff --git a/ode-evolve-benchmark.lua b/benchmarks/odeiv/odeiv-benchmark.lua index 0a11c098..8ed69cc0 100644 --- a/ode-evolve-benchmark.lua +++ b/benchmarks/odeiv/odeiv-benchmark.lua @@ -2,7 +2,7 @@ local template = require 'template' local ode_spec = {N = 2, eps_abs = 1e-6, eps_rel = 0, a_y = 1, a_dydt = 0} -local ode = template.load('rkf45.lua.in', ode_spec) +local ode = template.load('num/rkf45.lua.in', ode_spec) function f_vanderpol_gen(mu) return function(t, x, y) return y, -x + mu * y * (1-x^2) end diff --git a/ode-evolve-benchmark-vec.lua b/benchmarks/odeiv/odeiv-vec-benchmark.lua index 266366dc..022bb3da 100644 --- a/ode-evolve-benchmark-vec.lua +++ b/benchmarks/odeiv/odeiv-vec-benchmark.lua @@ -4,7 +4,7 @@ local template = require 'template' local ffi = require "ffi" local ode_spec = {N = 2, eps_abs = 1e-6, eps_rel = 0, a_y = 1, a_dydt = 0} -local ode = template.load('rkf45vec.lua.in', ode_spec) +local ode = template.load('num/rkf45vec.lua.in', ode_spec) function f_vanderpol_gen(mu) return function(t, y, f) diff --git a/examples/ode-jit.lua b/examples/ode-jit.lua new file mode 100644 index 00000000..9ceafa24 --- /dev/null +++ b/examples/ode-jit.lua @@ -0,0 +1,129 @@ + + -- ODE examples, examples/ode-jit.lua + -- + -- Copyright (C) 2009 Francesco Abbate + -- + -- This program is free software; you can redistribute it and/or modify + -- it under the terms of the GNU General Public License as published by + -- the Free Software Foundation; either version 3 of the License, or (at + -- your option) any later version. + -- + -- This program is distributed in the hope that it will be useful, but + -- WITHOUT ANY WARRANTY; without even the implied warranty of + -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + -- General Public License for more details. + -- + -- You should have received a copy of the GNU General Public License + -- along with this program; if not, write to the Free Software + -- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + -- + +local echo = gsl.echo +local sin, cos, exp = math.sin, math.cos, math.exp + +local template = require 'template' + +local ode_spec = {N = 2, eps_abs = 1e-6, eps_rel = 0, a_y = 1, a_dydt = 0} +local ode = template.load('num/rkf45.lua.in', ode_spec) + +function f_vanderpol_gen(mu) + return function(t, x, y) return y, -x + mu * y * (1-x^2) end +end + +local function xyodeplot(f, t0, t1, x0, y0, h0, tsmp) + local s = ode.new() + local evol = ode.evolve + local ln = graph.path(x0, y0) + ode.init(s, t0, h0, f, x0, y0) + local tsmp = 0.04 + for t = tsmp, t1, tsmp do + while s.t < t do + evol(s, f, t) + end + ln:line_to(s.y[0], s.y[1]) + end + + local p = graph.plot('ODE integration example') + p:addline(ln) + p:show() + return p +end + +local function modeplot(s, f, t0, y0, t1, tsmp) + local n = #y0 + local t = {} + for k=1, n do t[k] = graph.path(t0, y0[k]) end + local evol = ode.evolve + + if tsmp then + for t = tsmp, t1, tsmp do + while s.t < t do + evol(s, f, t) + end + for k=1, n do t[k]:line_to(s.t, s.y[k-1]) end + end + else + while s.t < t1 do + evol(s, f, t1) + for k=1, n do t[k]:line_to(s.t, s.y[k-1]) end + end + end + + local p = graph.plot('ODE integration example') + for k=1, n do p:addline(t[k], graph.rainbow(k)) end + p:show() + return p +end + +function demo1() + local odef = function(t, x, y) + return -y-x^2, 2*x - y^3 + end + + local x, y = 1, 1 + + return xyodeplot(odef, 0, 30, x, y, 1e-3, 0.04) +end + +function demo2() + local f = f_vanderpol_gen(10.0) + + local t0, t1, h = 0, 50, 0.01 + local y0 = {1, 0} + + local s = ode.new() + ode.init(s, t0, h, f, y0[1], y0[2]) + + return modeplot(s, f, t0, y0, t1) +end + +local function f_sincos_damp(damp_f) + return function (t, s, c) + return c - damp_f * s, -s - damp_f * c + end +end + +function demo3() + local damp_f = 0.04 + local x, y = 1, 0 + local s = ode.new() + local f = f_sincos_damp(damp_f) + local t0, t1, h0, hsamp = 0, 100, 1e-6, 0.5 + ode.init(s, t0, h0, f, x, y) + local ln = graph.path(t0, x) + for t= hsamp, t1, hsamp do + while s.t < t do + ode.evolve(s, f, t) + end + ln:line_to(s.t, s.y[0]) + end + local p = graph.plot() + p:addline(ln) + p:addline(graph.fxline(function(t) return cos(t) * exp(- damp_f * t) end, 0, 100), 'blue', {{'dash', 7,3}}) + p:show() + return p +end + +echo 'demo1() - ODE integration example' +echo 'demo2() - GSL example of Var der Pol oscillator integration' +echo 'demo3() - Examples of damped harmonic oscillator' diff --git a/examples/qag-bessel.lua b/examples/qag-bessel.lua index 360c01de..2ac47c30 100644 --- a/examples/qag-bessel.lua +++ b/examples/qag-bessel.lua @@ -1,10 +1,11 @@ local template = require 'template' -local q = template.load('num/qag.lua.in', {limit=64, order=21}) +local qag = template.load('num/qag.lua.in', {limit=64, order=21}) +local qng = template.load('num/qng.lua.in', {}) local sin, cos, pi = math.sin, math.cos, math.pi local epsabs, epsrel = 1e-6, 0.01 -function bessel_gen(n) +function bessel_gen(n, q) local xs local fint = function(t) return cos(n*t - xs*sin(t)) end return function(x) @@ -13,9 +14,17 @@ function bessel_gen(n) end end -local J4 = bessel_gen(4) +local J4 = bessel_gen(4, qag) +local J4b = bessel_gen(4, qng) -p = graph.fxplot(J4, 0, 30*pi) +w = graph.window 'v..' + +p1 = graph.plot('J4 bessel function') +p2 = graph.plot('J4 bessel function') +p1:addline(graph.fxline(J4, 0, 30*pi), 'red') +p2:addline(graph.fxline(J4b, 0, 8*pi), 'blue') +w:attach(p1, 1) +w:attach(p2, 2) local xold, xsmp = 0, 1 for x=0, 12*pi, 0.001 do diff --git a/examples/qng-bessel.lua b/examples/qng-bessel.lua deleted file mode 100644 index 8fbe1648..00000000 --- a/examples/qng-bessel.lua +++ /dev/null @@ -1,27 +0,0 @@ -local template = require 'template' -local qng = template.require 'num/qng' - -local sin, cos, pi = math.sin, math.cos, math.pi -local epsabs, epsrel = 1e-6, 1e-6 - -function bessel_gen(n) - local xs - local fint = function(t) return cos(n*t - xs*sin(t)) end - return function(x) - xs = x - return qng(fint, 0, pi, epsabs, epsrel) / pi - end -end - -local J4 = bessel_gen(4) - --- p = graph.fxplot(J4, 0, 30) - -local xold, xsmp = 0, 1 -for x=0, 12*pi, 0.001 do - local y = J4(x) - if x - xold > xsmp and y > -1 then - print(y) - xold = x - end -end diff --git a/num/gauss-kronrod-x-wgs.lua.in b/num/gauss-kronrod-x-wgs.lua.in index e712f9f3..b452af92 100644 --- a/num/gauss-kronrod-x-wgs.lua.in +++ b/num/gauss-kronrod-x-wgs.lua.in @@ -1,3 +1,4 @@ + # -- Gauss quadrature weights and kronrod quadrature abscissae and # -- weights as evaluated with 80 decimal digit arithmetic by # -- L. W. Fullerton, Bell Labs, Nov. 1981. diff --git a/ode-defs.lua.in b/num/ode-defs.lua.in index f94ba27f..f94ba27f 100644 --- a/ode-defs.lua.in +++ b/num/ode-defs.lua.in diff --git a/num/qag.lua.in b/num/qag.lua.in index 6c92c202..7dc80a2d 100644 --- a/num/qag.lua.in +++ b/num/qag.lua.in @@ -1,4 +1,25 @@ +# -- num/qag.lua.in +# -- +# -- Copyright (C) 2009-2011 Francesco Abbate +# -- +# -- This program is free software; you can redistribute it and/or modify +# -- it under the terms of the GNU General Public License as published by +# -- the Free Software Foundation; either version 3 of the License, or (at +# -- your option) any later version. +# -- +# -- This program is distributed in the hope that it will be useful, but +# -- WITHOUT ANY WARRANTY; without even the implied warranty of +# -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# -- General Public License for more details. +# -- +# -- You should have received a copy of the GNU General Public License +# -- along with this program; if not, write to the Free Software +# -- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +# -- + +# -- Adapted from the GSL Library, version 1.14 + # -- template parameters are 'limit' and 'order' # half_order = (order - 1) / 2 diff --git a/num/qng.lua b/num/qng.lua index 4c29de5e..46210567 100644 --- a/num/qng.lua +++ b/num/qng.lua @@ -1,8 +1,8 @@ --- integration/qng.lua +-- num/qng.lua -- -- Adapted to Lua by Francesco Abbate -- --- Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough +-- Copyright (C) 2009-2011 Francesco Abbate -- -- This program is free software; you can redistribute it and/or modify -- it under the terms of the GNU General Public License as published by diff --git a/num/qng.lua.in b/num/qng.lua.in index ea0aa23d..308eb177 100644 --- a/num/qng.lua.in +++ b/num/qng.lua.in @@ -1,6 +1,6 @@ -# -- integration/qng.c +# -- num/qng.lua.in # -- -# -- Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough +# -- Copyright (C) 2009-2011 Francesco Abbate # -- # -- This program is free software; you can redistribute it and/or modify # -- it under the terms of the GNU General Public License as published by diff --git a/rk4.lua.in b/num/rk4.lua.in index a06cec0e..1c82c787 100644 --- a/rk4.lua.in +++ b/num/rk4.lua.in @@ -1,7 +1,7 @@ # order = 4
-$(include 'ode-defs.lua.in')
+$(include 'num/ode-defs.lua.in')
local eps_abs, eps_rel, a_y, a_dydt = 1e-12, 1e-4, 1, 0
diff --git a/rkf45.lua.in b/num/rkf45.lua.in index 7e5ba926..d4012f00 100644 --- a/rkf45.lua.in +++ b/num/rkf45.lua.in @@ -1,9 +1,36 @@ +# -- num/rkf45.lua.in +# -- +# -- Copyright (C) 2009-2011 Francesco Abbate +# -- +# -- This program is free software; you can redistribute it and/or modify +# -- it under the terms of the GNU General Public License as published by +# -- the Free Software Foundation; either version 3 of the License, or (at +# -- your option) any later version. +# -- +# -- This program is distributed in the hope that it will be useful, but +# -- WITHOUT ANY WARRANTY; without even the implied warranty of +# -- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# -- General Public License for more details. +# -- +# -- You should have received a copy of the GNU General Public License +# -- along with this program; if not, write to the Free Software +# -- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +# -- + +# -- Adapted from the GSL Library, version 1.14 + +# -- Runge-Kutta-Fehlberg 4(5) + +# -- Reference eg. Hairer, E., Norsett S.P., Wanner, G. Solving ordinary +# -- differential equations I, Nonstiff Problems, 2nd revised edition, +# -- Springer, 2000. + local abs, max, min = math.abs, math.max, math.min # order = 5 -$(include 'ode-defs.lua.in') +$(include 'num/ode-defs.lua.in') # ah = { 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0 } # b3 = { 3.0/32.0, 9.0/32.0 } @@ -17,8 +44,8 @@ $(include 'ode-defs.lua.in') # c5 = -1371249.0/7618050.0 # c6 = 277020.0/7618050.0 --- These are the differences of fifth and fourth order coefficients --- for error estimation */ +# -- These are the differences of fifth and fourth order coefficients +# -- for error estimation # ec = { 0.0, 1.0 / 360.0, 0.0, -128.0 / 4275.0, -2197.0 / 75240.0, 1.0 / 50.0, 2.0 / 55.0 } diff --git a/rkf45vec.lua.in b/num/rkf45vec.lua.in index a477342b..a477342b 100644 --- a/rkf45vec.lua.in +++ b/num/rkf45vec.lua.in diff --git a/ode-test-vec.lua b/ode-test-vec.lua deleted file mode 100644 index f0348273..00000000 --- a/ode-test-vec.lua +++ /dev/null @@ -1,104 +0,0 @@ - -local gsl, graph = gsl or _G, graph or _G - -local template = require 'template' - -local ode_spec = {N = 2, eps_abs = 1e-6, eps_rel = 0, a_y = 1, a_dydt = 0} -local ode = template.load('rkf45vec.lua.in', ode_spec) - -local ffi = require "ffi" - -local sin, cos, exp, pi = math.sin, math.cos, math.exp, math.pi -local path, plot, fxplot, fxline = graph.path, graph.plot, graph.fxplot, graph.fxline - -function f_ode1(t, y, f) - f[0] = -y[1] - y[0]^2 - f[1] = 2*y[0] - y[1]^3 -end - -function f_vanderpol_gen(mu) - return function(t, y, f) - f[0] = y[1] - f[1] = -y[0] + mu * y[1] * (1-y[0]^2) - end -end - -local damp_f = 0.04 -function f_sincos_damp(t, y, f) - f[0] = y[1] - damp_f * y[0] - f[1] = -y[0] - damp_f * y[1] -end - -function test1() - local y = ffi.new('double[2]', {1, 0}) - local s = ode.new() - local f = f_sincos_damp - local t0, t1, h0, hsamp = 0, 100, 1e-6, 0.5 - ode.init(s, t0, h0, f, y) - local ln = path(t0, y[0]) - for t= hsamp, t1, hsamp do - while s.t < t do - ode.evolve(s, f, t) - end - ln:line_to(s.t, s.y[0]) - end - local p = plot() - p:addline(ln) - p:addline(fxline(function(t) return cos(t) * exp(- damp_f * t) end, 0, 100), 'blue', {{'dash', 7,3}}) - p:show() - return p -end - -function test2() - local s, f = ode.new(), f_ode1 - local y = ffi.new('double[2]', {1, 1}) - local t0, t1, h0, hsamp = 0, 30, 0.04, 0.04 - local ln = path(y[0], y[1]) - ode.init(s, t0, h0, f, y) - for t= hsamp, t1, hsamp do - while s.t < t do - ode.evolve(s, f, t) - end - ln:line_to(s.y[0], s.y[1]) - end - local p = plot() - p:addline(ln) - p:show() - return p -end - -function test3() - local f = f_vanderpol_gen(10.0) - local s = ode.new() - local y = ffi.new('double[2]', {1, 0}) - local t0, t1, h0 = 0, 200, 0.01 - for k=1, 10 do - ode.init(s, t0, h0, f, y) - while s.t < t1 do - ode.evolve(s, f, t1) - end - print(s.t, s.y[0], s.y[1]) - end -end - -function test4() - local f = f_vanderpol_gen(10.0) - local y = ffi.new('double[2]', {1, 0}) - local t0, t1, h0 = 0, 100, 1e-6 - - local s = ode.new() - ode.init(s, t0, h0, f, y) - - local lna, lnb = path(t0, y[0]), path(t0, y[1]) - while s.t < t1 do - ode.evolve(s, f, t1) - lna:line_to(s.t, s.y[0]) - lnb:line_to(s.t, s.y[1]) - end - - local p = plot() - p:addline(lna) - p:addline(lnb, 'blue') - p:show() - return p -end diff --git a/ode-test.lua b/ode-test.lua deleted file mode 100644 index be13665a..00000000 --- a/ode-test.lua +++ /dev/null @@ -1,97 +0,0 @@ - -local gsl, graph = gsl or _G, graph or _G - -local template = require 'template' - -local ode_spec = {N = 2, eps_abs = 1e-6, eps_rel = 0, a_y = 1, a_dydt = 0} -local ode = template.load('rkf45.lua.in', ode_spec) - -local sin, cos, exp, pi = math.sin, math.cos, math.exp, math.pi -local path, plot, fxplot, fxline = graph.path, graph.plot, graph.fxplot, graph.fxline - -function f_ode1(t, p, q) - return -q - p^2, 2*p - q^3 -end - -function f_vanderpol_gen(mu) - return function(t, x, y) return y, -x + mu * y * (1-x^2) end -end - -local damp_f = 0.04 -function f_sincos_damp(t, s, c) - return c - damp_f * s, -s - damp_f * c -end - -function test1() - local x, y = 1, 0 - local s = ode.new() - local f = f_sincos_damp - local t0, t1, h0, hsamp = 0, 100, 1e-6, 0.5 - ode.init(s, t0, h0, f, x, y) - local ln = path(t0, x) - for t= hsamp, t1, hsamp do - while s.t < t do - ode.evolve(s, f, t) - end - ln:line_to(s.t, s.y[0]) - end - local p = plot() - p:addline(ln) - p:addline(fxline(function(t) return cos(t) * exp(- damp_f * t) end, 0, 100), 'blue', {{'dash', 7,3}}) - p:show() - return p -end - -function test2() - local s, f = ode.new(), f_ode1 - local x, y = 1, 1 - local t0, t1, h0, hsamp = 0, 30, 0.04, 0.04 - local ln = path(x, y) - ode.init(s, t0, h0, f, x, y) - for t= hsamp, t1, hsamp do - while s.t < t do - ode.evolve(s, f, t) - end - ln:line_to(s.y[0], s.y[1]) - end - local p = plot() - p:addline(ln) - p:show() - return p -end - -function test3() - local f = f_vanderpol_gen(10.0) - local s = ode.new() - local x, y = 1, 0 - local t0, t1, h0 = 0, 20000, 0.01 - for k=1, 10 do - ode.init(s, t0, h0, f, x, y) - while s.t < t1 do - ode.evolve(s, f, t1) - end - print(s.t, s.y[0], s.y[1]) - end -end - -function test4() - local f = f_vanderpol_gen(10.0) - local x, y = 1, 0 - local t0, t1, h0 = 0, 5, 1e-6 - - local s = ode.new() - ode.init(s, t0, h0, f, x, y) - - local lna, lnb = path(t0, x), path(t0, y) - while s.t < t1 do - ode.evolve(s, f, t1) - lna:line_to(s.t, s.y[0]) - lnb:line_to(s.t, s.y[1]) - end - - local p = plot() - p:addline(lna) - p:addline(lnb, 'blue') - p:show() - return p -end diff --git a/rk4-ref-evolve.c b/rk4-ref-evolve.c deleted file mode 100644 index 7e57d9dd..00000000 --- a/rk4-ref-evolve.c +++ /dev/null @@ -1,70 +0,0 @@ -#include <stdio.h> -#include <gsl/gsl_errno.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_odeiv.h> - -int -func (double t, const double y[], double f[], - void *params) -{ - double mu = *(double *)params; - f[0] = y[1]; - f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1); - return GSL_SUCCESS; -} - -int -jac (double t, const double y[], double *dfdy, - double dfdt[], void *params) -{ - double mu = *(double *)params; - gsl_matrix_view dfdy_mat - = gsl_matrix_view_array (dfdy, 2, 2); - gsl_matrix * m = &dfdy_mat.matrix; - gsl_matrix_set (m, 0, 0, 0.0); - gsl_matrix_set (m, 0, 1, 1.0); - gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0); - gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0)); - dfdt[0] = 0.0; - dfdt[1] = 0.0; - return GSL_SUCCESS; -} - -int -main (void) -{ - const gsl_odeiv_step_type * T - = gsl_odeiv_step_rk8pd; - - gsl_odeiv_step * s - = gsl_odeiv_step_alloc (T, 2); - gsl_odeiv_control * c - = gsl_odeiv_control_y_new (1e-6, 0.0); - gsl_odeiv_evolve * e - = gsl_odeiv_evolve_alloc (2); - - double mu = 10; - gsl_odeiv_system sys = {func, jac, 2, &mu}; - - double t = 0.0, t1 = 100.0; - double h = 1e-6; - double y[2] = { 1.0, 0.0 }; - - while (t < t1) - { - int status = gsl_odeiv_evolve_apply (e, c, s, - &sys, - &t, t1, - &h, y); - - if (status != GSL_SUCCESS) - break; - - printf ("%.5e %.5e %.5e\n", t, y[0], y[1]); - } - - gsl_odeiv_evolve_free (e); - gsl_odeiv_control_free (c); - gsl_odeiv_step_free (s); - return 0; -} diff --git a/rk4-ref.c b/rk4-ref.c deleted file mode 100644 index 76e40275..00000000 --- a/rk4-ref.c +++ /dev/null @@ -1,60 +0,0 @@ - -#include <stdio.h> -#include <math.h> -#include <gsl/gsl_errno.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_odeiv.h> - -#define ODE_SYS_DIM 2 - -int -f_ode1 (double t, const double y[], double f[], void *params) -{ - double p = y[0], q = y[1]; - f[0] = - q - p*p; - f[1] = 2*p - q*q*q; - return GSL_SUCCESS; -} - -void -do_rk(gsl_odeiv_system *sys, double p0, double q0, double sample) -{ - gsl_odeiv_step *s = gsl_odeiv_step_alloc (gsl_odeiv_step_rk4, ODE_SYS_DIM); - double y[ODE_SYS_DIM], dydt[ODE_SYS_DIM], yerr[ODE_SYS_DIM]; - double t = 0, t1 = 2000, tsamp = 0, h0 = 0.001; - - y[0] = p0; - y[1] = q0; - - gsl_odeiv_step_apply (s, t, h0, y, yerr, NULL, dydt, sys); - t = t + h0; - - while (t < t1) - { - gsl_odeiv_step_apply (s, t, h0, y, yerr, dydt, dydt, sys); - t = t + h0; - if (sample > 0.0 && t - tsamp > sample) - { - printf("%g %g %g\n", t, y[0], y[1]); - tsamp = t; - } - } - - printf(">> %g %g %g\n", t, y[0], y[1]); - - gsl_odeiv_step_free (s); -} - -int -main() -{ - gsl_odeiv_system sys[1] = {{f_ode1, NULL, ODE_SYS_DIM, NULL}}; - int k; - for (k = 0; k < 10; k++) - { - double th = 3.14159265359 / 4; - double p0 = cos(th), q0 = sin(th); - do_rk (sys, p0, q0, -1.0); - } - return 0; -} diff --git a/rk4.lua b/rk4.lua deleted file mode 100644 index 93c392cf..00000000 --- a/rk4.lua +++ /dev/null @@ -1,118 +0,0 @@ - -ffi = require 'ffi' -darray = ffi.typeof("double[?]") - -ffi.cdef[[ - typedef struct { - double y[2]; - double dydt[2]; - } rk4_state; -]] - -function rk4_new() - return ffi.new('rk4_state') -end - -function rk4_step(y_0, y_1, dydt, h, t, f) - -- Makes a Runge-Kutta 4th order advance with step size h. - - -- initial values of variables y. - local y0_0, y0_1 = y_0, y_1 - - -- work space - local ytmp_0, ytmp_1 - - -- Runge-Kutta coefficients. Contains values of coefficient k1 - -- in the beginning - local k_0, k_1 = dydt[0], dydt[1] - - -- k1 step - y_0 = y_0 + h / 6 * k_0 - ytmp_0 = y0_0 + 0.5 * h * k_0 - y_1 = y_1 + h / 6 * k_1 - ytmp_1 = y0_1 + 0.5 * h * k_1 - - -- k2 step - k_0, k_1 = f(t + 0.5 * h, ytmp_0, ytmp_1) - - y_0 = y_0 + h / 3 * k_0 - ytmp_0 = y0_0 + 0.5 * h * k_0 - y_1 = y_1 + h / 3 * k_1 - ytmp_1 = y0_1 + 0.5 * h * k_1 - - -- k3 step - k_0, k_1 = f(t + 0.5 * h, ytmp_0, ytmp_1) - - y_0 = y_0 + h / 3 * k_0 - ytmp_0 = y0_0 + h * k_0 - y_1 = y_1 + h / 3 * k_1 - ytmp_1 = y0_1 + h * k_1 - - -- k4 step - k_0, k_1 = f(t + h, ytmp_0, ytmp_1) - - return y_0 + h / 6 * k_0, y_1 + h / 6 * k_1 -end - -function rk4_apply(s, t, h, yerr, f) - local y_0, y_1 = s.y[0], s.y[1] - local dydt = s.dydt - - -- First traverse h with one step (save to yonestep) - local yonestep_0, yonestep_1 = rk4_step (y_0, y_1, dydt, h, t, f) - - -- first step of h/2 - y_0, y_1 = rk4_step(y_0, y_1, dydt, h/2, t, f) - - dydt[0], dydt[1] = f(t + h/2, y_0, y_1) - - -- second step of h/2 - y_0, y_1 = rk4_step(y_0, y_1, dydt, h/2, t + h/2, f) - - -- Derivatives at output - dydt[0], dydt[1] = f(t + h, y_0, y_1) - - -- Error estimation - -- - -- yerr = C * 0.5 * | y(onestep) - y(twosteps) | / (2^order - 1) - -- - -- constant C is approximately 8.0 to ensure 90% of samples lie within - -- the error (assuming a gaussian distribution with prior p(sigma)=1/sigma.) - - yerr[0] = 4 * (y_0 - yonestep_0) / 15 - yerr[1] = 4 * (y_1 - yonestep_1) / 15 - - s.y[0], s.y[1] = y_0, y_1 -end - -function f_ode1(t, p, q) - return -q - p^2, 2*p - q^3 -end - -t0, t1, h0 = 0, 200, 0.001 - -function do_rk(p0, q0, sample) - local f = f_ode1 - local s = rk4_new() - local yerr = darray(2) - - s.y[0], s.y[1] = p0, q0 - s.dydt[0], s.dydt[1] = f(t, p0, q0) - - local t, tsamp = t0, t0 - while t < t1 do - rk4_apply(s, t, h0, yerr, f) - t = t + h0 - if sample and t - tsamp > sample then - print(t, s.y[0], s.y[1]) - tsamp = t - end - end - print(t, s.y[0], s.y[1]) -end - -for k=1, 10 do - local th = pi/4 -- *(k-1)/5 - local p0, q0 = cos(th), sin(th) - do_rk(p0, q0) -end diff --git a/rkf45vec-idamax.lua.in b/rkf45vec-idamax.lua.in deleted file mode 100644 index a02b6481..00000000 --- a/rkf45vec-idamax.lua.in +++ /dev/null @@ -1,181 +0,0 @@ - -local abs, max, min = math.abs, math.max, math.min - -# order = 5 - -# ah = { 1.0/4.0, 3.0/8.0, 12.0/13.0, 1.0, 1.0/2.0 } -# b3 = { 3.0/32.0, 9.0/32.0 } -# b4 = { 1932.0/2197.0, -7200.0/2197.0, 7296.0/2197.0} -# b5 = { 8341.0/4104.0, -32832.0/4104.0, 29440.0/4104.0, -845.0/4104.0} -# b6 = { -6080.0/20520.0, 41040.0/20520.0, -28352.0/20520.0, 9295.0/20520.0, -5643.0/20520.0} - -# c1 = 902880.0/7618050.0 -# c3 = 3953664.0/7618050.0 -# c4 = 3855735.0/7618050.0 -# c5 = -1371249.0/7618050.0 -# c6 = 277020.0/7618050.0 - -# -- These are the differences of fifth and fourth order coefficients -# -- for error estimation - -# ec = { 0.0, 1.0 / 360.0, 0.0, -128.0 / 4275.0, -2197.0 / 75240.0, 1.0 / 50.0, 2.0 / 55.0 } - -# y_err_only = (a_dydt == 0) - -local ffi = require "ffi" - -local vecsize = $(N) * ffi.sizeof('double') - -local function ode_new() - return ffi.new('ode_state') -end - -local function ode_init(s, t0, h0, f, y) - ffi.copy(s.y, y, vecsize) - f(t0, s.y, s.dydt) - s.t = t0 - s.h = h0 -end - -local function hadjust(rmax, h) - local S = 0.9 - if rmax > 1.1 then - local r = S / rmax^(1/$(order)) - r = max(0.2, r) - return r * h, -1 - elseif rmax < 0.5 then - local r = S / rmax^(1/($(order)+1)) - r = max(1, min(r, 5)) - return r * h, 1 - end - return h, 0 -end - -ffi.cdef[[ - typedef struct { - double t; - double h; - double y[$(N)]; - double dydt[$(N)]; - } ode_state; - - typedef struct { - double y0[$(N)]; - double ytmp[$(N)]; - double k1[$(N)]; - double k2[$(N)]; - double k3[$(N)]; - double k4[$(N)]; - double k5[$(N)]; - double k6[$(N)]; - double yerr[$(N)]; - } ode_workspace; - - void cblas_daxpy (const int N, const double ALPHA, - const double * X, const int INCX, - double * Y, const int INCY); - - int cblas_idamax (const int N, const double * X, const int INCX); - - void cblas_dscal (const int N, const double ALPHA, double * X, const int INCX); -]] - -local ws = ffi.new('ode_workspace') - -local function rkf45_evolve(s, f, t1) - local t, h = s.t, s.h - local hadj, inc - - ffi.copy (ws.y0, s.y, vecsize) - ffi.copy (ws.k1, s.dydt, vecsize) - - if t + h > t1 then h = t1 - t end - - while h > 0 do - local rmax - - do - ffi.copy (ws.ytmp, s.y, vecsize) - ffi.C.cblas_daxpy ($(N), h * $(ah[1]), ws.k1, 1, ws.ytmp, 1) - - -- k2 step - f(t + $(ah[1]) * h, ws.ytmp, ws.k2) - - ffi.copy (ws.ytmp, s.y, vecsize) - ffi.C.cblas_daxpy ($(N), h * $(b3[1]), ws.k1, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b3[2]), ws.k2, 1, ws.ytmp, 1) - - -- k3 step - f(t + $(ah[2]) * h, ws.ytmp, ws.k3) - - ffi.copy (ws.ytmp, s.y, vecsize) - ffi.C.cblas_daxpy ($(N), h * $(b4[1]), ws.k1, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b4[2]), ws.k2, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b4[3]), ws.k3, 1, ws.ytmp, 1) - - -- k4 step - f(t + $(ah[3]) * h, ws.ytmp, ws.k4) - - ffi.copy (ws.ytmp, s.y, vecsize) - ffi.C.cblas_daxpy ($(N), h * $(b5[1]), ws.k1, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b5[2]), ws.k2, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b5[3]), ws.k3, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b5[4]), ws.k4, 1, ws.ytmp, 1) - - -- k5 step - f(t + $(ah[4]) * h, ws.ytmp, ws.k5) - - ffi.copy (ws.ytmp, s.y, vecsize) - ffi.C.cblas_daxpy ($(N), h * $(b6[1]), ws.k1, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b6[2]), ws.k2, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b6[3]), ws.k3, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b6[4]), ws.k4, 1, ws.ytmp, 1) - ffi.C.cblas_daxpy ($(N), h * $(b6[5]), ws.k5, 1, ws.ytmp, 1) - - -- k6 step and final sum - -- since k2 is no more used we could use k2 to store k6 - f(t + $(ah[5]) * h, ws.ytmp, ws.k6) - - ffi.C.cblas_daxpy ($(N), h * $(c1), ws.k1, 1, s.y, 1) - ffi.C.cblas_daxpy ($(N), h * $(c3), ws.k3, 1, s.y, 1) - ffi.C.cblas_daxpy ($(N), h * $(c4), ws.k4, 1, s.y, 1) - ffi.C.cblas_daxpy ($(N), h * $(c5), ws.k5, 1, s.y, 1) - ffi.C.cblas_daxpy ($(N), h * $(c6), ws.k6, 1, s.y, 1) - -# if not y_err_only then - f(t + h, s.y, s.dydt) -# end - - - ffi.fill(ws.yerr, vecsize) - ffi.C.cblas_daxpy ($(N), h * $(ec[2]), ws.k1, 1, ws.yerr, 1) - ffi.C.cblas_daxpy ($(N), h * $(ec[4]), ws.k3, 1, ws.yerr, 1) - ffi.C.cblas_daxpy ($(N), h * $(ec[5]), ws.k4, 1, ws.yerr, 1) - ffi.C.cblas_daxpy ($(N), h * $(ec[6]), ws.k5, 1, ws.yerr, 1) - ffi.C.cblas_daxpy ($(N), h * $(ec[7]), ws.k6, 1, ws.yerr, 1) - - -- TODO TODO TODO: implement the complete algorithm - d0 = $(eps_abs) - ffi.C.cblas_dscal ($(N), 1/abs(d0), ws.yerr, 1) - - local imax = ffi.C.cblas_idamax ($(N), ws.yerr, 1) - rmax = abs(ws.yerr[imax]) - end - - hadj, inc = hadjust(rmax, h) - if inc >= 0 then break end - - ffi.copy(s.y, ws.y0, vecsize) - h = hadj - end - -# if y_err_only then - f(t + h, s.y, s.dydt) -# end - s.t = t + h - s.h = hadj - - return h -end - -return {new= ode_new, init= ode_init, evolve= rkf45_evolve} |