From octave-sources-request at bevo dot che dot wisc dot edu  Thu Oct  7 13:41:43 1999
Subject: ode integrator suite - v1.0 --> the source
From: Marc Compere <compere at mail dot utexas dot edu>
To: octave-sources at bevo dot che dot wisc dot edu
Date: Thu, 07 Oct 1999 13:40:55 -0500


________________Begin ode23.m________________
function [tout, yout] = ode23(F, tspan, y0, tol, trace, count)
% ode23 (v1.0) Integrates a system of ordinary differential equations
using
% 2nd order Runge-Kutta formulas with a 3rd order method that is
Simpson's 1/3 rule.
% This requires 3 function evaluations per integration step.
%
% The error estimate formula and slopes are from
% Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle,
McGraw-Hill, 1985
%
% Usage:
%         [tout, yout] = ode23(F, tspan, y0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative COLUMN-vector; yprime(i) =
dy(i)/dt.
% tspan - [ tstart, tfinal ]
% y0    - Initial value COLUMN-vector.
% tol   - The desired accuracy. (Default: tol = 1.e-3).
% trace - If nonzero, each step is printed. (Default: trace = 0).
% count - if nonzero, variable 'counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'counter' is incremented in here, not in the state-dot file
%         simply make 'counter' global in the file that calls ode23
%
% OUTPUT:
% tout  - Returned integration time points (column-vector).
% yout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 6 October 1999

if nargin < 6, count = 0; end
if nargin < 5, trace = 0; end
if nargin < 4, tol = 1.e-3; end

pow = 1/8;

% The 2(3) coefficients:
 a(1,1)=0;
 a(2,1)=1/2;
 a(3,1)=-1; a(3,2)=2;
 % 2nd order b-coefficients
 b2(1)=0; b2(2)=1;
 % 5th order b-coefficients
 b3(1)=1/6; b3(2)=2/3; b3(3)=1/6;
 for i=1:3
  c(i)=sum(a(i,:));
 end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmax = (tfinal - t)/2.5;
hmin = (tfinal - t)/1e12;
h = (tfinal - t)/200; % initial guess at a step size
y = y0(:);            % this always creates a column vector, y
tout = t;             % first output time
yout = y.';           % first output solution

if count==1,
 global counter
 if ~exist('counter'),counter=0; end
end % if count

if trace
 clc, t, h, y
end

% The main loop
   while (t < tfinal) & (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % compute the slopes
      k(:,1)=feval(F,t,y);
      k(:,2)=feval(F,t+c(2)*h,y+h*(a(2,1)*k(:,1)));
      k(:,3)=feval(F,t+c(3)*h,y+h*(a(3,1)*k(:,1)+a(3,2)*k(:,2)));

      % increment the counter
      if count==1,
       counter = counter + 3;
      end % if

      % compute the 2nd order estimate
      y2=y + h*b2(2)*k(:,2);
      % compute the 3rd order estimate
      y3=y + h*(b3(1)*k(:,1) + b3(2)*k(:,2) + b3(3)*k(:,3));

      % estimate the local truncation error
      gamma1 = y3 - y2;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');
      tau = tol*max(norm(y,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         y = y3;
         tout = [tout; t];
         yout = [yout; y.'];
      end
      if trace
         home, t, h, y
      end

      % Update the step size
      if delta ~= 0.0
         h = min(hmax, 0.8*h*(tau/delta)^pow);
      end
   end;

   if (t < tfinal)
      disp('Step size grew too small.')
      t, h, y
   end
________________End ode23.m________________

________________Begin ode45.m________________
function [tout, yout] = ode45(F, tspan, y0, tol, trace, count)
% ode45 (v1.0) integrates a system of ordinary differential equations
using
% 4th & 5th order embedded Runge-Kutta-Fehlberg formulas.
% This requires 6 function evaluations per integration step.
% This particular implementation uses the 5th order estimate for yout.
%
% The Fehlberg 4(5) coefficients are from a tableu in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential
Equations
% and Differential-Agebraic Equations, Society for Industrial and
Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% The error estimate formula and slopes are from
% Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle,
McGraw-Hill, 1985
%
% Usage:
%         [tout, yout] = ode45(F, tspan, y0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative COLUMN-vector; yprime(i) =
dy(i)/dt.
% tspan - [ tstart, tfinal ]
% y0    - Initial value COLUMN-vector.
% tol   - The desired accuracy. (Default: tol = 1.e-6).
% trace - If nonzero, each step is printed. (Default: trace = 0).
% count - if nonzero, variable 'counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'counter' is incremented in here, not in the state-dot file
%         simply make 'counter' global in the file that calls ode45
%
% OUTPUT:
% tout  - Returned integration time points (column-vector).
% yout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 6 October 1999

if nargin < 6, count = 0; end
if nargin < 5, trace = 0; end
if nargin < 4, tol = 1.e-6; end

pow = 1/8;

% The Fehlberg 4(5) coefficients:
 a_(1,1)=0;
 a_(2,1)=1/4;
 a_(3,1)=3/32; a_(3,2)=9/32;
 a_(4,1)=1932/2197; a_(4,2)=-7200/2197; a_(4,3)=7296/2197;
 a_(5,1)=439/216; a_(5,2)=-8; a_(5,3)=3680/513; a_(5,4)=-845/4104; 
 a_(6,1)=-8/27; a_(6,2)=2; a_(6,3)=-3544/2565; a_(6,4)=1859/4104;
a_(6,5)=-11/40; 
 % 4th order b-coefficients
 b4_(1)=25/216; b4_(2)=0; b4_(3)=1408/2565; b4_(4)=2197/4104;
b4_(5)=-1/5;
 % 5th order b-coefficients
 b5_(1)=16/135; b5_(2)=0; b5_(3)=6656/12825; b5_(4)=28561/56430;
b5_(5)=-9/50; b5_(6)=2/55;
 for i=1:6
  c_(i)=sum(a_(i,:));
 end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmax = (tfinal - t)/2.5;
hmin = (tfinal - t)/1e9;
h = (tfinal - t)/100; % initial guess at a step size
y = y0(:);            % this always creates a column vector, y
tout = t;             % first output time
yout = y.';           % first output solution

if count==1,
 global counter
 if ~exist('counter'),counter=0; end
end % if count

if trace
 clc, t, h, y
end

% The main loop
   while (t < tfinal) & (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % compute the slopes
      k_(:,1)=feval(F,t,y);
      k_(:,2)=feval(F,t+c_(2)*h,y+h*(a_(2,1)*k_(:,1)));
      k_(:,3)=feval(F,t+c_(3)*h,y+h*(a_(3,1)*k_(:,1)+a_(3,2)*k_(:,2)));
     
k_(:,4)=feval(F,t+c_(4)*h,y+h*(a_(4,1)*k_(:,1)+a_(4,2)*k_(:,2)+a_(4,3)*k_(:,3)));
     
k_(:,5)=feval(F,t+c_(5)*h,y+h*(a_(5,1)*k_(:,1)+a_(5,2)*k_(:,2)+a_(5,3)*k_(:,3)+a_(5,4)*k_(:,4)));
     
k_(:,6)=feval(F,t+c_(6)*h,y+h*(a_(6,1)*k_(:,1)+a_(6,2)*k_(:,2)+a_(6,3)*k_(:,3)+a_(6,4)*k_(:,4)+a_(6,5)*k_(:,5)));

      % increment the counter
      if count==1,
       counter = counter + 6;
      end % if

      % compute the 4th order estimate
      y4=y + h*(b4_(1)*k_(:,1) + b4_(3)*k_(:,3) + b4_(4)*k_(:,4) +
b4_(5)*k_(:,5));
      % compute the 5th order estimate
      y5=y + h*(b5_(1)*k_(:,1) + b5_(3)*k_(:,3) + b5_(4)*k_(:,4) +
b5_(5)*k_(:,5) + b5_(6)*k_(:,6));

      % estimate the local truncation error
      gamma1 = y5 - y4;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');       % actual error
      tau = tol*max(norm(y,'inf'),1.0); % allowable error

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         y = y5;
         tout = [tout; t];
         yout = [yout; y.'];
      end
      if trace
         home, t, h, y
      end

      % Update the step size
      if delta ~= 0.0
         h = min(hmax, 0.8*h*(tau/delta)^pow);
      end
   end;

   if (t < tfinal)
      disp('Step size grew too small.')
      t, h, y
   end
________________End ode45.m________________

________________Begin ode78.m________________
function [tout, yout] = ode78(F, tspan, y0, tol, trace, count)
% ode78 (v1.0) Integrates a system of ordinary differential equations
using
% 7th order formulas.
% This requires 13 function evaluations per integration step.
%
% More may be found in the original author's text containing numerous
% applications on ordinary and partial differential equations using
Matlab:
%
%     Howard Wilson and Louis Turcotte, 'Advanced Mathematics and 
%     Mechanics Applications Using MATLAB', 2nd Ed, CRC Press, 1997
%
%
% [tout, yout] = ode78(F, tspan, y0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative COLUMN-vector; yprime(i) =
dy(i)/dt.
% tspan - [ tstart, tfinal ]
% y0    - Initial value COLUMN-vector.
% tol   - The desired accuracy. (Default: tol = 1.e-6).
% trace - If nonzero, each step is printed. (Default: trace = 0).
% count - if nonzero, variable 'counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'counter' is incremented in here, not in the state-dot file
%         simply make 'counter' global in the file that calls ode78
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% yout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).

%   Daljeet Singh & Howard Wilson
%   Dept. Of Electrical Engg., The University of Alabama.
%   11-24-1988.
%
%   modified by:
%   Marc Compere
%   compere at mail dot utexas dot edu
%   6 October 1999


% The Fehlberg coefficients:
alpha_ = [ 2./27. 1/9 1/6 5/12 .5 5/6 1/6 2/3 1/3 1 0 1 ]';
beta_ = [ [  2/27  0  0   0   0  0  0  0  0  0  0   0  0  ]
[  1/36 1/12  0  0  0  0  0  0   0  0  0  0  0  ]
[  1/24  0  1/8  0  0  0  0  0  0  0  0  0  0 ]
[  5/12  0  -25/16  25/16  0  0  0  0  0  0   0  0  0  ]
[ .05   0  0  .25  .2  0  0  0  0  0  0  0  0 ]
[ -25/108  0  0  125/108  -65/27  125/54  0  0  0  0  0  0   0  ]
[ 31/300  0  0  0  61/225  -2/9  13/900  0  0  0   0  0  0  ]
[ 2  0  0  -53/6  704/45  -107/9  67/90  3  0  0  0  0  0  ]
[ -91/108  0  0  23/108  -976/135  311/54  -19/60  17/6  -1/12  0  0  0 
0 ]
[2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0
0 0]
[ 3/205  0   0  0   0    -6/41  -3/205   -3/41     3/41   6/41   0   0 
0 ]
[-1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 ...
51/82 33/164 12/41 0 1 0]...
]';
 chi_ = [ 0 0 0 0 0 34/105 9/35 9/35 9/280 9/280 0 41/840 41/840]';
 psi_ = [1  0  0  0  0  0  0  0  0  0  1 -1  -1 ]';
pow = 1/8;
if nargin < 6, count = 0; end
if nargin < 5, trace = 0; end
if nargin < 4, tol = 1.e-6; end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
% the following step parameters are used in ODE45
% hmax = (tfinal - t)/5;
% hmin = (tfinal - t)/20000;
% h = (tfinal - t)/100;
% The following parameters were taken because the integrator has
% higher order than ODE45. This choice is somewhat subjective.
hmax = (tfinal - t)/2.5;
%hmin = (tfinal - t)/10000;
hmin = (tfinal - t)/1000000000;
h = (tfinal - t)/50;
y = y0(:);
f = y*zeros(1,13);
tout = t;
yout = y.';
tau = tol * max(norm(y, 'inf'), 1);

if count==1,
 global counter
 if ~exist('counter'),counter=0;,end
end % if count

if trace
%  clc, t, h, y
   clc, t, y
end
% The main loop
   while (t < tfinal) & (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes
      f(:,1) = feval(F,t,y);
      for j = 1: 12
         f(:,j+1) = feval(F, t+alpha_(j)*h, y+h*f*beta_(:,j));
      end

      % increment the counter
      if count==1,
       counter = counter + 13;
      end % if

      % Truncation error term
      gamma1 = h*41/840*f*psi_;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');
      tau = tol*max(norm(y,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         y = y + h*f*chi_;
         tout = [tout; t];
         yout = [yout; y.'];
      end
      if trace
         home, t, h, y
%        home, t, y
      end

      % Update the step size
      if delta ~= 0.0
         h = min(hmax, 0.8*h*(tau/delta)^pow);
      end
   end;

   if (t < tfinal)
      disp('SINGULARITY LIKELY.')
      t
   end
________________End ode78.m________________

________________Begin rk2fixed.m________________
function [tout, yout] = rk2fixed(F, tspan, y0, N, trace, count)
% rk2fixed (v1.0) integrates a system of ordinary differential equations
using a
% 2nd order Runge-Kutta formula called Ralston's method.
% This choice of 2nd order coefficients provides a minimum bound on
truncation error.
% For more, see Ralston & Rabinowitz (1978) or 
% Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle,
McGraw-Hill, 1985
%
% rk2fixed() requires 2 function evaluations per integration step.
%
% Usage:
%         [tout, yout] = rk2fixed(F, tspan, y0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem derivatives.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative COLUMN-vector; yprime(i) =
dy(i)/dt.
% tspan - [ tstart, tfinal ]
% y0    - Initial value COLUMN-vector.
% N     - number of steps used to span [ tstart, tfinal ]
% trace - If nonzero, each step is printed. (Default: trace = 0).
% count - if nonzero, variable 'counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'counter' is incremented in here, not in the state-dot file
%         simply make 'counter' global in the file that calls rk2fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% yout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 6 October 1999

if nargin < 6, count = 0; end
if nargin < 5, trace = 0; end
if nargin < 4, tol = 1.e-6; end

if count==1,
 global counter
 if ~exist('counter'),counter=0;,end
end % if count

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
hmax = (tfinal - t)/2.5;
hmin = (tfinal - t)/1e12;
h = (tfinal - t)/200; % initial guess at a step size
yout(1,:) = y0';      % this always creates a row vector, y, and is the
first output
tout(1)   = t0;       % first output time
t = t0;
y = y0(:);


if trace
 clc, t, h, y
end

% The main loop
h = (tspan(2)-tspan(1))/N;
neqs = size(y0(:),1);

for i=2:N+1,
     k1 = feval(F,t,y);
     k2 = feval(F,t+3/4*h,y+3/4*h*k1);

     % increment the counter
     if count==1,
      counter = counter + 2;
     end % if

     yout(i,:) = (y+h*(1/3*k1+2/3*k2))';
     tout(i) = t+h;
     y = yout(i,:)';
     t = tout(i);

     if trace,
      home, t, h, yout(i,:)'
     end

end
________________End rk2fixed.m________________

________________Begin rk4fixed.m________________
function [xout,yout] = rk4fixed(fxy,xspan,y0,N,trace,count)

% rk4fixed (v1.0) is a numerical integration routine.
% It requires 4 function evaluations per step.
%
% Usage:
%         [tout, yout] = rk4fixed(F, tspan, y0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem derivatives.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative COLUMN-vector; yprime(i) =
dy(i)/dt.
% tspan - [ tstart, tfinal ]
% y0    - Initial value COLUMN-vector.
% N     - number of steps used to span [ tstart, tfinal ]
% trace - if nonzero, each step is printed. (Default: trace = 0).
% count - if nonzero, variable 'counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'counter' is incremented in here, not in the state-dot file
%         simply make 'counter' global in the file that calls rk4fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% yout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 6 October 1999

if nargin < 5, count = 0; end
if nargin < 4, trace = 0; end

if count==1,
 global counter
 if ~exist('counter'),counter=0;,end
end % if count

x0 = xspan(1);
h = (xspan(2)-xspan(1))/N;
halfh = 0.5*h;
neqs = size(y0(:),1);
yout = zeros(N,neqs(1));
xout = zeros(1,N);
yout(1,:) = y0';
xout(1) = x0;
x = x0;
y = y0(:);

for i=2:N+1,
     RK1 = feval(fxy,x,y);
     xhalf = x+halfh;
     ytemp = y+halfh*RK1;
     RK2 = feval(fxy,xhalf,ytemp);
     ytemp = y+halfh*RK2;
     RK3 = feval(fxy,xhalf,ytemp);
     xfull = x+h;
     ytemp = y+h*RK3;
     RK4 = feval(fxy,xfull,ytemp);

     % increment the counter
     if count==1,
      counter = counter + 4;
     end % if

     yout(i,:) = (y+h/6*(RK1+2.0*(RK2+RK3)+RK4))';
     xout(i) = xfull;
     y = yout(i,:)';
     x = xout(i);

     if trace,
      home, x, h, y
     end

end
________________End rk4fixed.m________________

________________Begin rk8fixed.m________________
function [tout,yout] = rk8fixed(fxy,tspan,y0,N,trace,count)

% rk8fixed (v1.0) is a numerical integration routine.
% It requires 13 function evaluations per step.
%
% Usage:
%         [tout, yout] = rk8fixed(F, tspan, y0, N, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem derivatives.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative COLUMN-vector; yprime(i) =
dy(i)/dt.
% tspan - [ tstart, tfinal ]
% y0    - Initial value COLUMN-vector.
% N     - number of steps used to span [ tstart, tfinal ]
% trace - if nonzero, each step is printed. (Default: trace = 0).
% count - if nonzero, variable 'counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'counter' is incremented in here, not in the state-dot file
%         simply make 'counter' global in the file that calls rk4fixed
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% yout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 6 October 1999

if nargin < 5, count = 0; end
if nargin < 4, trace = 0; end

if count==1,
 global counter
 if ~exist('counter'),counter=0;,end
end % if count

alpha_ = [ 2./27. 1/9 1/6 5/12 .5 5/6 1/6 2/3 1/3 1 0 1 ]';
beta_ = [ [  2/27  0  0   0   0  0  0  0  0  0  0   0  0  ]
[  1/36 1/12  0  0  0  0  0  0   0  0  0  0  0  ]
[  1/24  0  1/8  0  0  0  0  0  0  0  0  0  0 ]
[  5/12  0  -25/16  25/16  0  0  0  0  0  0   0  0  0  ]
[ .05   0  0  .25  .2  0  0  0  0  0  0  0  0 ]
[ -25/108  0  0  125/108  -65/27  125/54  0  0  0  0  0  0   0  ]
[ 31/300  0  0  0  61/225  -2/9  13/900  0  0  0   0  0  0  ]
[ 2  0  0  -53/6  704/45  -107/9  67/90  3  0  0  0  0  0  ]
[ -91/108  0  0  23/108  -976/135  311/54  -19/60  17/6  -1/12  0  0  0 
0 ]
[2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0
0 0]
[ 3/205  0   0  0   0    -6/41  -3/205   -3/41     3/41   6/41   0   0 
0 ]
[-1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 ...
51/82 33/164 12/41 0 1 0]...
]';
chi_ = [ 0 0 0 0 0 34/105 9/35 9/35 9/280 9/280 0 41/840 41/840]';

t0 = tspan(1);
h = (tspan(2)-tspan(1))/N;
neqs = size(y0(:),1);
yout(1,:) = y0';
tout(1) = t0;
t = t0;
y = y0(:);
f = y*zeros(1,13);

for i=2:N+1,

     % Compute the slopes
     f(:,1) = feval(fxy,t,y);
     for j = 1:12
        f(:,j+1) = feval(fxy, t+alpha_(j)*h, y+h*f*beta_(:,j));
     end

     % increment the counter
     if count==1,
      counter = counter + 13;
     end % if

     t = t + h;
     y = y + h*f*chi_;
     tout = [tout; t];
     yout = [yout; y.'];

     if trace,
      home, t, h, y
     end

end
________________End rk8fixed.m________________

________________Begin pendulum.m________________
% This integrates a simple pendulum example with 6 different
integrators.
%
% As of now, ode78.m, ode45.m, ode23.m, rk8fixed.m, rk4fixed.m, and
rk2fixed.m.
% all produce output column vectors similar to Matlab.
%
% All integrators work in octave 2.0.14 and Matlab 5.2 with no
modification.
%
% Marc Compere
% compere at mail dot utexas dot edu
% 6 October 1999

clear

% allow global access to the parameters m, g, l, & b:
global m g l b counter

m=1;    % (kg)
g=9.81; % (m/s^2)
l=1.;   % (m)
b=2.;   % ((N-s)/m))

% integrator setup:
trace = 0;   % this is a (1/0) flag that puts output to the screen or
not
count = 1;   % this is a (1/0) flag that causes rk4fixed to increment
'counter' or not
counter = 0; % 'counter' is the number of right-hand-side function
evaluations, z_dot=f(z)
t0=0;
tfinal = 5; % (seconds)
tspan = [t0,tfinal];
% Initial Conditions: theta(t=0)=30(deg) & initially at rest
IC = [ 30*pi/180 0]'; % (rad), (rad/s)
sps = 100;            % sps -> step per second
Nsteps=(tfinal-t0)*sps;    % this creates sps number of integration
steps per second
tolerance = 1e-3;



t_begin_calcs=cputime;

   disp('Integrating using rk2fixed...')
   [t1,zrk2fixed] = rk2fixed('penddot',tspan,IC,Nsteps,trace,count);    
% fixed step integration
   rk2_counter = counter
   counter=0;

t(1)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using rk4fixed...')
   [t2,zrk4fixed] = rk4fixed('penddot',tspan,IC,Nsteps,trace,count);    
% fixed step integration
   rk4_counter = counter
   counter=0;

t(2)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using rk8fixed...')
   [t3,zrk8fixed] = rk8fixed('penddot',tspan,IC,Nsteps,trace,count);    
% fixed step integration
   rk8_counter = counter
   counter=0;

t(3)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode23...')
   [t4,zode23] = ode23('penddot',tspan,IC,tolerance,trace,count); % rk45
variable step integration
   ode23_counter = counter
   counter=0;

t(4)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode45...')
   [t5,zode45] = ode45('penddot',tspan,IC,tolerance,trace,count); % rk45
variable step integration
   ode45_counter = counter
   counter=0;

t(5)=cputime-t_begin_calcs;
t_begin_calcs=cputime;

   disp('Integrating using ode78...')
   [t6,zode78] = ode78('penddot',tspan,IC,tolerance,trace,count); % rk78
variable step integration
   ode78_counter = counter
   counter=0;

t(6)=cputime-t_begin_calcs;



disp('Elapsed times for each solver to integrate a pendulum:')
t



% plot that baby
%figure(1)
if ishold~=1, hold, end
plot(t1,zrk2fixed)
plot(t2,zrk4fixed)
plot(t3,zrk8fixed)
plot(t4,zode23)
plot(t5,zode45)
plot(t6,zode78)
if ishold==1, hold, end
________________End pendulum.m________________

________________Begin penddot.m________________
function zdot=penddot(t,z)
%
% This is an example derivative function file made to work
% with octave or Matlab.
% The physical system is a simple pendulum with damping.
%
% create is with the form shown below, then use it
% like this:
%    [t,z] = ode45('penddot',to,tfinal,IC,tolerance); 
%
% z is the state column vector and meant to be used only within this
m-file
% This function is meant to return the derivatives of the state variable
% given the state vector and time.
%
% Structure:	zdot = [ z1dot, z2dot, ... zNdot ]'
%
% eg.  ml^2*thetadd + b*thetad + m*g*l*sin(theta) = 0
%
%	z(1) = theta
%	z(2) = thetad 	( = z(1)dot )
%
% Convention: the lowest order states are first columnwise

global m g l b counter index

zdot=[ z(2) , 1/(m*l^2)*(-b*z(2)-m*g*l*sin(z(1)))]';

% remember to return a column vector
________________End penddot.m________________



-- 
_________________________________________________
Marc Compere, The University of Texas at Austin
e-mail: compere at mail dot utexas dot edu
www   : http://nerdlab.me.utexas.edu/~compere
office: ETC 4.134a, (512) 471-7347
_________________________________________________