| author | francesco <francesco.bbt@gmail.com> | 2011年03月15日 12:55:20 +0100 |
|---|---|---|
| committer | francesco <francesco.bbt@gmail.com> | 2011年03月15日 12:55:20 +0100 |
| commit | 60582ef229f9a57d7f7fd107cd9894a689200a11 (patch) | |
| tree | 160d144048b0032e2e35d929a75a42f90fe661da /matrix.lua | |
| parent | b2d6e7aae94c9d4e4f581594d5de6fd43158c81f (diff) | |
| download | gsl-shell-60582ef229f9a57d7f7fd107cd9894a689200a11.tar.gz | |
| -rw-r--r-- | matrix.lua | 71 |
diff --git a/matrix.lua b/matrix.lua index 8d7e6319..d4f4e08b 100644 --- a/matrix.lua +++ b/matrix.lua @@ -235,74 +235,3 @@ function gsl.linfit(gener, x, y, w) end return f, c end - -local function nlinfitwrapf(fmodel, x, y) - local n = gsl.dim(y) - return function(p, f, J) - for k=1, n do - local ym = fmodel(p, x[k], J and J:row(k)) - if f then f:set(k, 1, ym - y[k]) end - end - end - -end - -function gsl.nlinfit(f, x, y, p0) - local N = gsl.dim(y) - local P = gsl.dim(p0) - local s = gsl.nlfsolver {fdf= nlinfitwrapf(f, x, y), n= N, p0= p0} - s:run() - return function(x) return f(s.p, x) end, s.p -end - -function gsl.ode(spec) - local required = {N= 'number', eps_abs= 'number'} - local defaults = {eps_rel = 0, a_y = 1, a_dydt = 0} - local is_known = {rkf45= true, rk8pd= true} - - for k, tp in pairs(required) do - if type(spec[k]) ~= tp then - error(string.format('parameter %s should be a %s', k, tp)) - end - end - for k, v in pairs(defaults) do - if not spec[k] then spec[k] = v end - end - - local method = spec.method and spec.method or 'rkf45' - if not is_known[method] then error('unknown ode method: ' .. method) end - spec.method = nil - - local ode = template.load(string.format('num/%s.lua.in', method), spec) - - local ode_methods = { - evolve = function(s, f, t) - s._sync = false - return ode.evolve(s._state, f, t) - end, - init = function(s, t0, h0, f, ...) - s._sync = false - return ode.init(s._state, t0, h0, f, ...) - end - } - - local ODE = {__index= function(s, k) - if k == 't' then return s._state.t end - if k == 'y' then - if not s._sync then - for k=1, s.dim do - s._y[k] = s._state.y[k-1] - end - s._sync = true - end - return s._y - end - return ode_methods[k] - end} - - local dim = spec.N - local solver = {_state = ode.new(), dim= dim, _y = gsl.new(dim,1)} - setmetatable(solver, ODE) - - return solver -end |