|
0
|
1 function [tout,xout] = rk8fixed(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 % rk8fixed.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 % rk8fixed.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 % rk8fixed (v1.07) is an 8th order Runge-Kutta numerical integration routine. |
|
|
18 % It requires 13 function evaluations per step. This is not the most |
|
|
19 % efficient 8th order implementation. It was just the easiest to put |
|
|
20 % together as a variant from ode78.m. |
|
|
21 % |
|
|
22 % Usage: |
|
|
23 % [tout, xout] = rk8fixed(F, tspan, x0, Nsteps, ode_fcn_format, trace, count) |
|
|
24 % |
|
|
25 % INPUT: |
|
|
26 % F - String containing name of user-supplied problem derivatives. |
|
|
27 % Call: xprime = fun(t,x) where F = 'fun'. |
|
|
28 % t - Time or independent variable (scalar). |
|
|
29 % x - Solution column-vector. |
|
|
30 % xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt. |
|
|
31 % tspan - [ tstart, tfinal ] |
|
|
32 % x0 - Initial value COLUMN-vector. |
|
|
33 % Nsteps - number of steps used to span [ tstart, tfinal ] |
|
|
34 % ode_fcn_format - this specifies if the user-defined ode function is in |
|
|
35 % the form: xprime = fun(t,x) (ode_fcn_format=0, default) |
|
|
36 % or: xprime = fun(x,t) (ode_fcn_format=1) |
|
|
37 % Matlab's solvers comply with ode_fcn_format=0 while |
|
|
38 % Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1. |
|
|
39 % trace - If nonzero, each step is printed. (optional, default: trace = 0). |
|
|
40 % count - if nonzero, variable 'rhs_counter' is initalized, made global |
|
|
41 % and counts the number of state-dot function evaluations |
|
|
42 % 'rhs_counter' is incremented in here, not in the state-dot file |
|
|
43 % simply make 'rhs_counter' global in the file that calls rk4fixed |
|
|
44 % |
|
|
45 % OUTPUT: |
|
|
46 % tout - Returned integration time points (row-vector). |
|
|
47 % xout - Returned solution, one solution column-vector per tout-value. |
|
|
48 % |
|
|
49 % The result can be displayed by: plot(tout, xout). |
|
|
50 % |
|
|
51 % Marc Compere |
|
|
52 % compere@mail.utexas.edu |
|
|
53 % created : 06 October 1999 |
|
|
54 % modified: 15 May 2000 |
|
|
55 |
|
|
56 if nargin < 7, count = 0; end |
|
|
57 if nargin < 6, trace = 0; end |
|
|
58 if nargin < 5, Nsteps = 50/(tspan(2)-tspan(1)); end % <-- 50 is a guess for a default, |
|
|
59 % try verifying the solution with ode78 |
|
|
60 if nargin < 4, ode_fcn_format = 0; end |
|
|
61 |
|
|
62 if count==1, |
|
|
63 global rhs_counter |
|
|
64 if ~exist('rhs_counter'),rhs_counter=0;,end |
|
|
65 end % if count |
|
|
66 |
|
|
67 alpha_ = [ 2./27. 1/9 1/6 5/12 .5 5/6 1/6 2/3 1/3 1 0 1 ]'; |
|
|
68 beta_ = [ [ 2/27 0 0 0 0 0 0 0 0 0 0 0 0 ] |
|
|
69 [ 1/36 1/12 0 0 0 0 0 0 0 0 0 0 0 ] |
|
|
70 [ 1/24 0 1/8 0 0 0 0 0 0 0 0 0 0 ] |
|
|
71 [ 5/12 0 -25/16 25/16 0 0 0 0 0 0 0 0 0 ] |
|
|
72 [ .05 0 0 .25 .2 0 0 0 0 0 0 0 0 ] |
|
|
73 [ -25/108 0 0 125/108 -65/27 125/54 0 0 0 0 0 0 0 ] |
|
|
74 [ 31/300 0 0 0 61/225 -2/9 13/900 0 0 0 0 0 0 ] |
|
|
75 [ 2 0 0 -53/6 704/45 -107/9 67/90 3 0 0 0 0 0 ] |
|
|
76 [ -91/108 0 0 23/108 -976/135 311/54 -19/60 17/6 -1/12 0 0 0 0 ] |
|
|
77 [2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0 0 0] |
|
|
78 [ 3/205 0 0 0 0 -6/41 -3/205 -3/41 3/41 6/41 0 0 0 ] |
|
|
79 [-1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 ... |
|
|
80 51/82 33/164 12/41 0 1 0]... |
|
|
81 ]'; |
|
|
82 chi_ = [ 0 0 0 0 0 34/105 9/35 9/35 9/280 9/280 0 41/840 41/840]'; |
|
|
83 |
|
|
84 % Initialization |
|
|
85 t = tspan(1); |
|
|
86 h = (tspan(2)-tspan(1))/Nsteps; |
|
|
87 xout(1,:) = x0'; |
|
|
88 tout(1) = t; |
|
|
89 x = x0(:); |
|
|
90 f = x*zeros(1,13); |
|
|
91 |
|
|
92 for i=1:Nsteps, |
|
|
93 |
|
|
94 % Compute the slopes |
|
|
95 if (ode_fcn_format==0), |
|
|
96 f(:,1) = feval(F,t,x); |
|
|
97 for j = 1:12 |
|
|
98 f(:,j+1) = feval(F, t+alpha_(j)*h, x+h*f*beta_(:,j)); |
|
|
99 end |
|
|
100 else, |
|
|
101 f(:,1) = feval(F,x,t); |
|
|
102 for j = 1:12 |
|
|
103 f(:,j+1) = feval(F, x+h*f*beta_(:,j), t+alpha_(j)*h); |
|
|
104 end |
|
|
105 end % if (ode_fcn_format==0) |
|
|
106 |
|
|
107 % increment rhs_counter |
|
|
108 if count==1, |
|
|
109 rhs_counter = rhs_counter + 13; |
|
|
110 end % if |
|
|
111 |
|
|
112 t = t + h; |
|
|
113 x = x + h*f*chi_; |
|
|
114 tout = [tout; t]; |
|
|
115 xout = [xout; x.']; |
|
|
116 |
|
|
117 if trace, |
|
|
118 home, t, h, x |
|
|
119 end |
|
|
120 |
|
|
121 end |
