Re: [OctDev] Ode23s

[davide p. <dav...@ho...>]

function [tout,xout] = ode23s(FUN,tspan,x0,options)

% This file is intended for use with Octave.
% ode23s.m 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 2, or (at your option)
% any later version.
%
% ode23.m 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 at www.gnu.org/copyleft/gpl.html.
%
% --------------------------------------------------------------------


if nargin < 4, options=odeset();end

% Initialization
d=1/(2+ sqrt(2));
a=1/2;
e32=6+sqrt(2);

if size(options.RelTol)~=size([]) %user-defined relative tolerance
tol=options.RelTol;
else
tol = 1.e-3;
end

t = tspan(1);
tfinal = tspan(2);

if size(options.MaxStep)~=size([]) %user-defined max step size
hmax=options.MaxStep;
else
hmax = (tfinal - t)/2.5;
end

hmin = 16*eps*(tfinal - t);

if size(options.InitialStep)~=size([]) %user-defined initial step size
h=options.InitialStep;
else
h = (tfinal - t)/200; % initial guess at a step size
end 

x = x0(:);            % this always creates a column vector, x
tout = t;             % first output time
xout = x.';           % first output solution


% The main loop
   while (t < tfinal) && (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end       
      
      % Jacobian matrix, dfxpdp
    
      J=dfxpdp(t,x,FUN);               
      T=(feval(FUN,x,t+hmin)-feval(FUN,x,t))/hmin;
         
      % Wolfbrandt coefficient

      W=eye(length(x0))-h*d*J;


      
      % compute the slopes

         F(:,1)=feval(FUN,x,t);
         k(:,1)=W\(F(:,1)+ h*d*T);
         F(:,2)=feval(FUN, x+a*h*k(:,1),t+a*h);
         k(:,2)=W\((F(:,2) - k(:,1))) + k(:,1);
         
         % compute the 2nd order estimate
         x2=x + h*k(:,2);
         
         % 3rd order, needed in error forumula
         F(:,3)=feval(FUN,x2,t+h);
         k(:,3)=W\(F(:,3)-e32*(k(:,2)-F(:,2))-2*(k(:,1)-F(:,1))+h*d*T);

      % estimate the error
      error = (h/6)*(norm(k(:,1)-2*k(:,2)+k(:,3)));

      % Estimate the acceptable error
      tau = tol*max(norm(x,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if error <= tau
         
         t = t + h;
         x = x2;  %no local extrapolation, FSAL
         tout = [tout; t];

         if size(options.Mass)~=size([]) %user-defined mass matrix

         M=options.Mass;
         xout = [xout;(M\x).'];  
     
         else

         xout = [xout; x.'];

         end

      % Update the step size
      
      if error == 0.0
       error = 1e-16;
      end

      h = min(hmax, h*1.25);    %auto-adaptive step update

      else
      
      h = max(hmin, h*0.5);    %auto-adaptive step update

      end
      


   end;
   
   if (t < tfinal)
      disp('Step size grew too small.')
      t, h, x
   end