Mercurial > forge

comparison extra/ode/ode45.m @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children a0d3391e59e2
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 function [tout, xout] = ode45(F,tspan,x0,ode_fcn_format,tol,trace,count)
2
3 % Copyright (C) 2000 Marc Compere
4 % This file is intended for use with Octave.
5 % ode45.m is free software; you can redistribute it and/or modify it
6 % under the terms of the GNU General Public License as published by
7 % the Free Software Foundation; either version 2, or (at your option)
8 % any later version.
9 %
10 % ode45.m is distributed in the hope that it will be useful, but
11 % WITHOUT ANY WARRANTY; without even the implied warranty of
12 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 % General Public License for more details at www.gnu.org/copyleft/gpl.html.
14 %
15 % --------------------------------------------------------------------
16 %
17 % ode45 (v1.07) integrates a system of ordinary differential equations using
18 % 4th & 5th order embedded Runge-Kutta-Fehlberg formulas.
19 % This requires 6 function evaluations per integration step.
20 % This particular implementation uses the 5th order estimate for xout, although
21 % the truncation error is of order(h^4), therefore this method, overall is
22 % a 4th-order method.
23 %
24 % The Fehlberg 4(5) coefficients are from a tableu in
25 % U.M. Ascher, L.R. Petzold, Computer Methods for Ordinary Differential Equations
26 % and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
27 % (SIAM), Philadelphia, 1998
28 %
29 % The error estimate formula and slopes are from
30 % Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle, McGraw-Hill, 1985
31 %
32 % Usage:
33 % [tout, xout] = ode45(F, tspan, x0, ode_fcn_format, tol, trace, count)
34 %
35 % INPUT:
36 % F - String containing name of user-supplied problem description.
37 % Call: xprime = fun(t,x) where F = 'fun'.
38 % t - Time (scalar).
39 % x - Solution column-vector.
40 % xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
41 % tspan - [ tstart, tfinal ]
42 % x0 - Initial value COLUMN-vector.
43 % ode_fcn_format - this specifies if the user-defined ode function is in
44 % the form: xprime = fun(t,x) (ode_fcn_format=0, default)
45 % or: xprime = fun(x,t) (ode_fcn_format=1)
46 % Matlab's solvers comply with ode_fcn_format=0 while
47 % Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
48 % tol - The desired accuracy. (optional, default: tol = 1.e-6).
49 % trace - If nonzero, each step is printed. (optional, default: trace = 0).
50 % count - if nonzero, variable 'rhs_counter' is initalized, made global
51 % and counts the number of state-dot function evaluations
52 % 'rhs_counter' is incremented in here, not in the state-dot file
53 % simply make 'rhs_counter' global in the file that calls ode45
54 %
55 % OUTPUT:
56 % tout - Returned integration time points (column-vector).
57 % xout - Returned solution, one solution column-vector per tout-value.
58 %
59 % The result can be displayed by: plot(tout, xout).
60 %
61 % Marc Compere
62 % compere@mail.utexas.edu
63 % created : 06 October 1999
64 % modified: 15 May 2000
65
66 if nargin < 7, count = 0; end
67 if nargin < 6, trace = 0; end
68 if nargin < 5, tol = 1.e-6; end
69 if nargin < 4, ode_fcn_format = 0; end
70
71 pow = 1/8;
72
73 % The Fehlberg 4(5) coefficients:
74 a_(1,1)=0;
75 a_(2,1)=1/4;
76 a_(3,1)=3/32; a_(3,2)=9/32;
77 a_(4,1)=1932/2197; a_(4,2)=-7200/2197; a_(4,3)=7296/2197;
78 a_(5,1)=439/216; a_(5,2)=-8; a_(5,3)=3680/513; a_(5,4)=-845/4104;
79 a_(6,1)=-8/27; a_(6,2)=2; a_(6,3)=-3544/2565; a_(6,4)=1859/4104; a_(6,5)=-11/40;
80 % 4th order b-coefficients
81 b4_(1)=25/216; b4_(2)=0; b4_(3)=1408/2565; b4_(4)=2197/4104; b4_(5)=-1/5;
82 % 5th order b-coefficients
83 b5_(1)=16/135; b5_(2)=0; b5_(3)=6656/12825; b5_(4)=28561/56430; b5_(5)=-9/50; b5_(6)=2/55;
84 for i=1:6
85 c_(i)=sum(a_(i,:));
86 end
87
88 % Initialization
89 t0 = tspan(1);
90 tfinal = tspan(2);
91 t = t0;
92 hmax = (tfinal - t)/2.5;
93 hmin = (tfinal - t)/1e9;
94 h = (tfinal - t)/100; % initial guess at a step size
95 x = x0(:); % this always creates a column vector, x
96 tout = t; % first output time
97 xout = x.'; % first output solution
98
99 if count==1,
100 global rhs_counter
101 if ~exist('rhs_counter'),rhs_counter=0; end
102 end % if count
103
104 if trace
105 clc, t, h, x
106 end
107
108 % The main loop
109 while (t < tfinal) & (h >= hmin)
110 if t + h > tfinal, h = tfinal - t; end
111
112 % compute the slopes
113 if (ode_fcn_format==0),
114 k_(:,1)=feval(F,t,x);
115 k_(:,2)=feval(F,t+c_(2)*h,x+h*(a_(2,1)*k_(:,1)));
116 k_(:,3)=feval(F,t+c_(3)*h,x+h*(a_(3,1)*k_(:,1)+a_(3,2)*k_(:,2)));
117 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)));
118 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)));
119 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)));
120 else,
121 k_(:,1)=feval(F,x,t);
122 k_(:,2)=feval(F,x+h*(a_(2,1)*k_(:,1)),t+c_(2)*h);
123 k_(:,3)=feval(F,x+h*(a_(3,1)*k_(:,1)+a_(3,2)*k_(:,2)),t+c_(3)*h);
124 k_(:,4)=feval(F,x+h*(a_(4,1)*k_(:,1)+a_(4,2)*k_(:,2)+a_(4,3)*k_(:,3)),t+c_(4)*h);
125 k_(:,5)=feval(F,x+h*(a_(5,1)*k_(:,1)+a_(5,2)*k_(:,2)+a_(5,3)*k_(:,3)+a_(5,4)*k_(:,4)),t+c_(5)*h);
126 k_(:,6)=feval(F,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)),t+c_(6)*h);
127 end % if (ode_fcn_format==0)
128
129 % increment rhs_counter
130 if count==1,
131 rhs_counter = rhs_counter + 6;
132 end % if
133
134 % compute the 4th order estimate
135 x4=x + h*(b4_(1)*k_(:,1) + b4_(3)*k_(:,3) + b4_(4)*k_(:,4) + b4_(5)*k_(:,5));
136 % compute the 5th order estimate
137 x5=x + h*(b5_(1)*k_(:,1) + b5_(3)*k_(:,3) + b5_(4)*k_(:,4) + b5_(5)*k_(:,5) + b5_(6)*k_(:,6));
138
139 % estimate the local truncation error
140 gamma1 = x5 - x4;
141
142 % Estimate the error and the acceptable error
143 delta = norm(gamma1,'inf'); % actual error
144 tau = tol*max(norm(x,'inf'),1.0); % allowable error
145
146 % Update the solution only if the error is acceptable
147 if delta <= tau
148 t = t + h;
149 x = x5; % <-- using the higher order estimate is called 'local extrapolation'
150 tout = [tout; t];
151 xout = [xout; x.'];
152 end
153 if trace
154 home, t, h, x
155 end
156
157 % Update the step size
158 if delta == 0.0
159 delta = 1e-16;
160 end
161 h = min(hmax, 0.8*h*(tau/delta)^pow);
162
163 end;
164
165 if (t < tfinal)
166 disp('Step size grew too small.')
167 t, h, x
168 end