From octave-sources-request at bevo dot che dot wisc dot edu  Thu Oct 21 14:35:43 1999
Subject: ode integrator suite - v1.01 --> 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, 21 Oct 1999 14:35:46 -0500

This is a multi-part message in MIME format.
--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

Octave users,

	Changes include the fix (or addition) of a default number of steps for
the fixed-step integrators.  Thanks to Rolf Fabian
<fabian at tu-cottbus dot de> for pointing out the bug.

	In addition, the whole set of 6 integrators has been cleaned up and
given more consistent variable names throughout.

Regards,
Marc
_________________________________________________
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
_________________________________________________
--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="ode23.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="ode23.m"

function [tout, xout] = ode23(F, tspan, x0, tol, trace, count)
% ode23 (v1.01) Integrates a system of ordinary differential equations using
% 2nd & 3rd order Runge-Kutta formulas.  The particular 3rd order method is
% Simpson's 1/3 rule.
% This particular implementation uses the 3rd order estimate for xout, although
% the truncation error is of order(h^2), therefore this method, overall is
% a 2nd-order method.
% 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, xout] = ode23(F, tspan, x0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where F = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - 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).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 21 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
x = x0(:);            % this always creates a column vector, x
tout = t;             % first output time
xout = x.';           % first output solution

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

if trace
 clc, t, h, x
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,x);
      k(:,2)=feval(F,t+c(2)*h,x+h*(a(2,1)*k(:,1)));
      k(:,3)=feval(F,t+c(3)*h,x+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
      x2=x + h*b2(2)*k(:,2);
      % compute the 3rd order estimate
      x3=x + h*(b3(1)*k(:,1) + b3(2)*k(:,2) + b3(3)*k(:,3));

      % estimate the local truncation error
      gamma1 = x3 - x2;

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

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         x = x3;    % <-- using the higher order estimate is called 'local extrapolation'
         tout = [tout; t];
         xout = [xout; x.'];
      end
      if trace
         home, t, h, x
      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, x
   end

--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="ode45.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="ode45.m"

function [tout, xout] = ode45(F, tspan, x0, 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 xout, although
% the truncation error is of order(h^4), therefore this method, overall is
% a 4th-order method.
%
% 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, xout] = ode45(F, tspan, x0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where F = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - 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).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 21 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
x = x0(:);            % this always creates a column vector, x
tout = t;             % first output time
xout = x.';           % first output solution

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

if trace
 clc, t, h, x
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,x);
      k_(:,2)=feval(F,t+c_(2)*h,x+h*(a_(2,1)*k_(:,1)));
      k_(:,3)=feval(F,t+c_(3)*h,x+h*(a_(3,1)*k_(:,1)+a_(3,2)*k_(:,2)));
      k_(:,4)=feval(F,t+c_(4)*h,x+h*(a_(4,1)*k_(:,1)+a_(4,2)*k_(:,2)+a_(4,3)*k_(:,3)));
      k_(:,5)=feval(F,t+c_(5)*h,x+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,x+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
      x4=x + h*(b4_(1)*k_(:,1) + b4_(3)*k_(:,3) + b4_(4)*k_(:,4) + b4_(5)*k_(:,5));
      % compute the 5th order estimate
      x5=x + 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 = x5 - x4;

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

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         x = x5;    % <-- using the higher order estimate is called 'local extrapolation'
         tout = [tout; t];
         xout = [xout; x.'];
      end
      if trace
         home, t, h, x
      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, x
   end

--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="ode78.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="ode78.m"

function [tout, xout] = ode78(F, tspan, x0, tol, trace, count)
% ode78 (v1.0) Integrates a system of ordinary differential equations using
% 7th order formulas.
% This particular implementation uses the 8th order estimate for xout, which
% is called 'local extrapolation'.  The truncation error, gamma1, is of order(h^8).
% Therefore this method, overall is a 7th-order method.
% 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, xout] = ode78(F, tspan, x0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where F = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - 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).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).

%   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
%   21 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;
x = x0(:);
f = x*zeros(1,13);
tout = t;
xout = x.';
tau = tol * max(norm(x, 'inf'), 1);

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

if trace
%  clc, t, h, x
   clc, t, x
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,x);
      for j = 1: 12
         f(:,j+1) = feval(F, t+alpha_(j)*h, x+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(x,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         x = x + h*f*chi_;  % this integrator uses local extrapolation
         tout = [tout; t];
         xout = [xout; x.'];
      end
      if trace
         home, t, h, x
%        home, t, x
      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

--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="penddot.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="penddot.m"

function xdot=penddot(t,x)
%
% 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,x] = ode45('penddot',to,tfinal,IC,tolerance); 
%
% x 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:	xdot = [ x1dot, x2dot, ... xNdot ]'
%
% eg.  ml^2*thetadd + b*thetad + m*g*l*sin(theta) = 0
%
%	x(1) = theta
%	x(2) = thetad 	( = x(1)dot )
%
% Convention: the lowest order states are first columnwise

global m g l b counter index

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

% remember to return a column vector

--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="pendulum.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="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
% 21 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, x_dot=f(x)
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,xrk2fixed] = 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,xrk4fixed] = 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,xrk8fixed] = 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,xode23] = 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,xode45] = 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,xode78] = 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,xrk2fixed)
plot(t2,xrk4fixed)
plot(t3,xrk8fixed)
plot(t4,xode23)
plot(t5,xode45)
plot(t6,xode78)
if ishold==1, hold, end




--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="readme.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="readme.txt"

This is the readme.txt file for octave_ode_solvers version 1.01.

This directory contains 8 files that provide an octave user with
several options for numerically integrating ode.
There are 3 fixed-step Runge-Kutta algorithms and
3 variable step Runge-Kutta-Fehlberg algorithms.
All are explicit RK formulas that work well with nonstiff problems.

----------------------------------------------------------------------
The archive octave_ode_solvers_v1.0.tar.gz should contain:

   - ode23.m    : variable step, 2nd-3rd order, Runge-Kutta, single-step method
   - ode45.m    : variable step, 4th-5th order, Runge-Kutta, single-step method
   - ode78.m    : variable step, 7th-8th order, Runge-Kutta, single-step method

   - rk2fixed.m : fixed step, 2nd order, Runge-Kutta, single-step method
   - rk4fixed.m : fixed step, 4th order, Runge-Kutta, single-step method
   - rk8fixed.m : fixed step, 8th order, Runge-Kutta, single-step method

   - pendulum.m : a sample m-file script that runs all solvers
   - penddot.m  : derivative function file, returning dy/dt for a simple pendulum
   - readme.txt : this file

----------------------------------------------------------------------
Steps for testing these ode solvers in octave version 2.0.14,
from a unix shell:

(1) unzip and untar the archive:
	gunzip octave_ode_solvers.tar.gz
	tar xvf octave_ode_solvers.tar
(2) change directories into the newly created directory and start octave
	cd octave_ode_solvers
	octave
(3) run the sample pemdulum script that will sequentially run all integrators
	pendulum

----------------------------------------------------------------------
I've made an effort to make these portable to most octave installations
as well as for use in most Matlab versions.  These work with no modification
in Matlab v5.2.  If you want to use these in Matlab, however,
you'll do yourself a favor by renaming ode45.m and ode23.m to something else,
like ode45_octave.m and ode23_octave.m.  This is because Matlab already has
two integrators named ode45 and ode23.
Don't forget to change the function names inside ode45_octave.m and
ode23_octave.m as well.

A little effort has been made to create ode45 and ode23 with similar
argument structures to their Matlab counterparts.  Only the most basic function
calls to ode45 and ode23 match up in both Matlab's integrators and these.
Feel free to change them as you see fit.  You are welcome to mail me for useful
changes I (may) make and continue to distribute from the same ftp site.

In general, the higher the integration order, the smaller the local 
error is at each time step.  This is demonstrated by ode78 generating far fewer
steps for solving the same problem over the same time interval with the same error
criterion as ode23.
The cost of the higher order solvers is the number of function evaluations.

   - ode23.m    : requires 3 function evaluations per step
   - ode45.m    : requires 6 function evaluations per step
   - ode78.m    : requires 13 function evaluations per step

   - rk2fixed.m : requires 2 function evaluations per step
   - rk4fixed.m : requires 4 function evaluations per step
   - rk8fixed.m : requires 13 function evaluations per step

If you get to looking at pendulum.m, try turning on the 'trace' variable
for screen output.  This slows things down, though.

----------------------------------------------------------------------
Two of the original files rk4fixed.m and ode78.m came from other
people.  rk4fixed.m was written by:
     Dr. Raul Longoria, Dept. of Mechanical Engineering,
     The Univ. of Texas at Austin and
ode78.m was originally written by:
     Dr. Howard Wilson & 
     Daljeet Singh, Dept. Of Electrical Engineering,
     The University of Alabama

The other files were created from the structure of these 
along with coefficients from standard numerical methods books.
Dr. Longoria has given permission to redistribute rk4fixed.m.

ode78.m was origianlly found at:
ftp://ftp.mathworks.com/pub/contrib/v4/diffeq/ode78.m.
& I am redistributing a modified version with Dr. Wilson's
permission.  

Numerous applications in ordinary and partial differential equations
can be found in Dr. Wilson's text:

     Howard Wilson and Louis Turcotte, 'Advanced Mathematics and 
     Mechanics Applications Using MATLAB', 2nd Ed, CRC Press, 1997

----------------------------------------------------------------------

Bugs or changes or comments should be directed to the email address
below.  The latest and greatest versions will hopefully befound
on the octave-source list at:
	http://www.che.wisc.edu/octave/mailing-lists/octave-sources/1999/
but should always be found at:
	http://marc.me.utexas.edu/tmp/octave_ode_solvers/


Marc Compere
compere at mail dot utexas dot edu
21 October 1999

--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="rk2fixed.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="rk2fixed.m"

function [tout, xout] = rk2fixed(F, tspan, x0, 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, xout] = rk2fixed(F, tspan, x0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem derivatives.
%         Call: xprime = fun(t,x) where F = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - 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).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 21 October 1999

if nargin < 6, count = 0; end
if nargin < 5, trace = 0; end
if nargin < 4, N = 400/(tspan(2)-tspan(1)); end % <-- 400 is a guess for a default,
                                                %  try verifying the solution with rk4fixed

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

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/N;
xout(1,:) = x0';           
tout(1)   = t;
x = x0(:);

if trace
 clc, t, h, x
end

% The main loop
h = (tspan(2)-tspan(1))/N;

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

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

     t = t + h;
     x = (x+h*(1/3*k1+2/3*k2));
     tout = [tout; t];
     xout = [xout; x.'];

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

end

--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="rk4fixed.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="rk4fixed.m"

function [tout,xout] = rk4fixed(F,tspan,x0,N,trace,count)

% rk4fixed (v1.01) is a numerical integration routine.
% It requires 4 function evaluations per step.
%
% Usage:
%         [tout, xout] = rk4fixed(F, tspan, x0, tol, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem derivatives.
%         Call: xprime = fun(t,x) where F = 'fun'.
%         t      - Time or independent variable (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - 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).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 21 October 1999

if nargin < 6, count = 0; end
if nargin < 5, trace = 0; end
if nargin < 4, N = 100/(tspan(2)-tspan(1)); end % <-- 100 is a guess for a default,
                                                %  try verifying the solution with rk8fixed

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

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/N;
xout(1,:) = x0';
tout(1) = t;
x = x0(:);
halfh = 0.5*h;

for i=2:N+1,
     RK1 = feval(F,t,x);
     thalf = t+halfh;
     xtemp = x+halfh*RK1;
     RK2 = feval(F,thalf,xtemp);
     xtemp = x+halfh*RK2;
     RK3 = feval(F,thalf,xtemp);
     tfull = t+h;
     xtemp = x+h*RK3;
     RK4 = feval(F,tfull,xtemp);

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

     t = t + h;
     x = (x+h/6*(RK1+2.0*(RK2+RK3)+RK4));
     tout = [tout; t];
     xout = [xout; x.'];

     if trace,
      home, t, h, x
     end

end

--------------9142EFF7F8699162710D3FAE
Content-Type: text/plain; charset=us-ascii;
 name="rk8fixed.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="rk8fixed.m"

function [tout,xout] = rk8fixed(F,tspan,x0,N,trace,count)

% rk8fixed (v1.01) is a numerical integration routine.
% It requires 13 function evaluations per step.  This is not the most
% efficient 8th order implementation.  It was just the easiest to put
% together as a variant from ode78.m.
%
% Usage:
%         [tout, xout] = rk8fixed(F, tspan, x0, N, trace, count)
%
% INPUT:
% F     - String containing name of user-supplied problem derivatives.
%         Call: xprime = fun(t,x) where F = 'fun'.
%         t      - Time or independent variable (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - 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).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).
%
% Marc Compere
% compere at mail dot utexas dot edu
% 21 October 1999

if nargin < 6, count = 0; end
if nargin < 5, trace = 0; end
if nargin < 4, N = 50/(tspan(2)-tspan(1)); end % <-- 50 is a guess for a default,
                                                %  try verifying the solution with ode78

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]';

% Initialization
t = tspan(1);
h = (tspan(2)-tspan(1))/N;
xout(1,:) = x0';
tout(1) = t;
x = x0(:);
f = x*zeros(1,13);

for i=2:N+1,

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

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

     t = t + h;
     x = x + h*f*chi_;
     tout = [tout; t];
     xout = [xout; x.'];

     if trace,
      home, t, h, x
     end

end

--------------9142EFF7F8699162710D3FAE--