Mercurial > forge

comparison extra/ode/rk2fixed.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] = rk2fixed(F,tspan,x0,Nsteps,ode_fcn_format,trace,count)
2
3 % Copyright (C) 2000 Marc Compere
4 % This file is intended for use with Octave.
5 % rk2fixed.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 % rk2fixed.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 % rk2fixed (v1.07) integrates a system of ordinary differential equations using a
18 % 2nd order Runge-Kutta formula called Ralston's method.
19 % This choice of 2nd order coefficients provides a minimum bound on truncation error.
20 % For more, see Ralston & Rabinowitz (1978) or
21 % Numerical Methods for Engineers, 2nd Ed., Chappra & Cannle, McGraw-Hill, 1985
22 %
23 % rk2fixed() requires 2 function evaluations per integration step.
24 %
25 % Usage:
26 % [tout, xout] = rk2fixed(F, tspan, x0, Nsteps, ode_fcn_format, trace, count)
27 %
28 % INPUT:
29 % F - String containing name of user-supplied problem derivatives.
30 % Call: xprime = fun(t,x) where F = 'fun'.
31 % t - Time (scalar).
32 % x - Solution column-vector.
33 % xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
34 % tspan - [ tstart, tfinal ]
35 % x0 - Initial value COLUMN-vector.
36 % Nsteps - number of steps used to span [ tstart, tfinal ]
37 % ode_fcn_format - this specifies if the user-defined ode function is in
38 % the form: xprime = fun(t,x) (ode_fcn_format=0, default)
39 % or: xprime = fun(x,t) (ode_fcn_format=1)
40 % Matlab's solvers comply with ode_fcn_format=0 while
41 % Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
42 % trace - If nonzero, each step is printed. (optional, default: trace = 0).
43 % count - if nonzero, variable 'rhs_counter' is initalized, made global
44 % and counts the number of state-dot function evaluations
45 % 'rhs_counter' is incremented in here, not in the state-dot file
46 % simply make 'rhs_counter' global in the file that calls rk2fixed
47 %
48 % OUTPUT:
49 % tout - Returned integration time points (row-vector).
50 % xout - Returned solution, one solution column-vector per tout-value.
51 %
52 % The result can be displayed by: plot(tout, xout).
53 %
54 % Marc Compere
55 % compere@mail.utexas.edu
56 % created : 06 October 1999
57 % modified: 15 May 2000
58
59 if nargin < 7, count = 0; end
60 if nargin < 6, trace = 0; end
61 if nargin < 5, Nsteps = 400/(tspan(2)-tspan(1)); end % <-- 400 is a guess for a default,
62 % try verifying the solution with rk4fixed
63 if nargin < 4, ode_fcn_format = 0; end
64
65 if count==1,
66 global rhs_counter
67 if ~exist('rhs_counter'),rhs_counter=0;,end
68 end % if count
69
70 % Initialization
71 t = tspan(1);
72 h = (tspan(2)-tspan(1))/Nsteps;
73 xout(1,:) = x0';
74 tout(1) = t;
75 x = x0(:);
76
77 if trace
78 clc, t, h, x
79 end
80
81 % The main loop
82 h = (tspan(2)-tspan(1))/Nsteps;
83
84 for i=1:Nsteps,
85 if (ode_fcn_format==0),
86 k1 = feval(F,t,x);
87 k2 = feval(F,t+3/4*h,x+3/4*h*k1);
88 else,
89 k1 = feval(F,x,t);
90 k2 = feval(F,x+3/4*h*k1,t+3/4*h);
91 end % if (ode_fcn_format==0)
92
93 % increment rhs_counter
94 if count==1,
95 rhs_counter = rhs_counter + 2;
96 end % if
97
98 t = t + h;
99 x = (x+h*(1/3*k1+2/3*k2));
100 tout = [tout; t];
101 xout = [xout; x.'];
102
103 if trace,
104 home, t, h, x
105 end
106
107 end