Mercurial > octave-nkf
diff scripts/ode/private/runge_kutta_45_dorpri.m @ 20635:a22d8a2eb0e5
fix adaptive strategy in ode solvers.
* script/ode/ode45.m: remove unused option OutputSave
* script/ode/private/integrate_adaptive.m: rewrite algorithm
to be more compatible.
* script/ode/private/runge_kutta_45_dorpri.m: use kahan summation
for time increment.
author | Carlo de Falco <carlo.defalco@polimi.it> |
---|---|
date | Sun, 11 Oct 2015 18:44:58 +0200 |
parents | b7ac1e94266e |
children |
line wrap: on
line diff
--- a/scripts/ode/private/runge_kutta_45_dorpri.m Sat Oct 10 16:52:59 2015 -0700 +++ b/scripts/ode/private/runge_kutta_45_dorpri.m Sun Oct 11 18:44:58 2015 +0200 @@ -58,7 +58,8 @@ ## @seealso{odepkg} function [t_out, x_out, x_est, k] = ... - runge_kutta_45_dorpri (f, t, x, dt, varargin) + runge_kutta_45_dorpri (f, t, x, dt, opts = [], k_vals = [], + t_out = t + dt) persistent a = [0 0 0 0 0 0; 1/5 0 0 0 0 0; @@ -80,17 +81,18 @@ aa = dt * a; k = zeros (rows (x), 7); - if (nargin >= 5) # options are passed - args = varargin{1}.vfunarguments; - if (nargin >= 6) # both the options and the k values are passed - k(:,1) = varargin{2}(:,end); # FSAL property - else - k(:,1) = feval (f, t, x, args{:}); - endif + if (! isempty (opts)) # extra arguments for function evaluator + args = opts.vfunarguments; else args = {}; endif + if (! isempty (k_vals)) # k values from previous step are passed + k(:,1) = k_vals(:,end); # FSAL property + else + k(:,1) = feval (f, t, x, args{:}); + endif + k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).', args{:}); k(:,3) = feval (f, s(3), x + k(:,1:2) * aa(3, 1:2).', args{:}); k(:,4) = feval (f, s(4), x + k(:,1:3) * aa(4, 1:3).', args{:}); @@ -98,13 +100,13 @@ k(:,6) = feval (f, s(6), x + k(:,1:5) * aa(6, 1:5).', args{:}); ## compute new time and new values for the unknowns - t_out = t + dt; + ## t_out = t + dt; x_out = x + k(:,1:6) * cc(:); # 5th order approximation ## if the estimation of the error is required if (nargout >= 3) ## new solution to be compared with the previous one - k(:,7) = feval (f, t + dt, x_out, args{:}); + k(:,7) = feval (f, t_out, x_out, args{:}); cc_prime = dt * c_prime; x_est = x + k * cc_prime(:); endif