Mercurial > octave-nkf

comparison scripts/ode/private/runge_kutta_45_dorpri.m @ 20572:72cd24aa5f7a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up and vetorize Dormant&Prince RK timestepper * scripts/ode/private/runge_kutta_45_dorpri: cleanup and vectorize.
author Carlo de Falco <carlo.defalco@polimi.it>
date Fri, 02 Oct 2015 17:48:58 +0200
parents 6256f6e366ac
children 3339c9bdfe6a
comparison
equal deleted inserted replaced
20571:d0f886a030b7 20572:72cd24aa5f7a
1 ## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> 1 ## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it>
2 ## Copyright (C) 2015, Carlo de Falco
2 ## 3 ##
3 ## This file is part of Octave. 4 ## This file is part of Octave.
4 ## 5 ##
5 ## Octave is free software; you can redistribute it and/or modify it 6 ## Octave 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 ## under the terms of the GNU General Public License as published by
62 ## 63 ##
63 ## @seealso{odepkg} 64 ## @seealso{odepkg}
64 65
65 function varargout = runge_kutta_45_dorpri (f, t, x, dt, varargin) 66 function varargout = runge_kutta_45_dorpri (f, t, x, dt, varargin)
66 67
67 options = varargin{1}; 68 persistent a = [0 0 0 0 0 0;
68 k = zeros (rows (x), 6); 69 1/5 0 0 0 0 0;
70 3/40 9/40 0 0 0 0;
71 44/45 -56/15 32/9 0 0 0;
72 19372/6561 -25360/2187 64448/6561 -212/729 0 0;
73 9017/3168 -355/33 46732/5247 49/176 -5103/18656 0;
74 35/384 0 500/1113 125/192 -2187/6784 11/84];
75 persistent b = [0 1/5 3/10 4/5 8/9 1 1];
76 persistent c = [(35/384) 0 (500/1113) (125/192) (-2187/6784) (11/84)];
77 ## x_est according to Shampine 1986:
78 ## persistent c_prime = [(1951/21600) 0 (22642/50085) (451/720), ...
79 ## (-12231/42400) (649/6300) (1/60)];
80 persistent c_prime = [(5179/57600) 0 (7571/16695) (393/640), ...
81 (-92097/339200) (187/2100) (1/40)];
82
83 s = t + dt * b;
84 cc = dt * c;
85 aa = dt * a;
86
87 args = varargin{1}.vfunarguments;
88 k = zeros (rows (x), 7);
69 89
70 if (nargin == 5) # only the options are passed 90 if (nargin == 5) # only the options are passed
71 k(:,1) = feval (f, t , x, options.vfunarguments{:}); 91 k(:,1) = feval (f, t, x, args{:});
72 elseif (nargin == 6) # both the options and the k values are passed 92 elseif (nargin == 6) # both the options and the k values are passed
73 k(:,1) = varargin{2}(:,end); # FSAL property 93 k(:,1) = varargin{2}(:,end); # FSAL property
74 endif 94 endif
75 95
76 k(:,1) = feval (f, t, x, options.vfunarguments{:}); 96 k(:,1) = feval (f, s(1), x, args{:});
77 k(:,2) = feval (f, t + (1/5)*dt, ... 97 k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).', args{:});
78 x + dt * (1/5)*k(:,1), ... 98 k(:,3) = feval (f, s(3), x + k(:,1:2) * aa(3, 1:2).', args{:});
79 options.vfunarguments{:}); 99 k(:,4) = feval (f, s(4), x + k(:,1:3) * aa(4, 1:3).', args{:});
80 k(:,3) = feval (f, t + (3/10)*dt, ... 100 k(:,5) = feval (f, s(5), x + k(:,1:4) * aa(5, 1:4).', args{:});
81 x + dt * ((3/40)*k(:,1) + (9/40)*k(:,2)), ... 101 k(:,6) = feval (f, s(6), x + k(:,1:5) * aa(6, 1:5).', args{:});
82 options.vfunarguments{:});
83 k(:,4) = feval (f, t + (4/5)*dt, ...
84 x + dt * ((44/45)*k(:,1) - (56/15)*k(:,2) + (32/9)*k(:,3)), ...
85 options.vfunarguments{:});
86 k(:,5) = feval (f, t + (8/9)*dt, ...
87 x + dt * ((19372/6561)*k(:,1) - (25360/2187)*k(:,2) ...
88 + (64448/6561)*k(:,3) - (212/729)*k(:,4)), ...
89 options.vfunarguments{:});
90 k(:,6) = feval (f, t + dt, ...
91 x + dt * ((9017/3168)*k(:,1) - (355/33)*k(:,2) ...
92 + (46732/5247)*k(:,3) + (49/176)*k(:,4) ...
93 - (5103/18656)*k(:,5)), ...
94 options.vfunarguments{:});
95 102
96 ## compute new time and new values for the unkwnowns 103 ## compute new time and new values for the unkwnowns
97 varargout{1} = t + dt; 104 varargout{1} = t + dt;
98 varargout{2} = x + dt * ((35/384)*k(:,1) + (500/1113)*k(:,3) 105 varargout{2} = x + k(:,1:6) * cc(:); # 5th order approximation
99 + (125/192)*k(:,4) - (2187/6784)*k(:,5)
100 + (11/84)*k(:,6)); # 5th order approximation
101 106
102 ## if the estimation of the error is required 107 ## if the estimation of the error is required
103 if (nargout >= 3) 108 if (nargout >= 3)
104 ## new solution to be compared with the previous one 109 ## new solution to be compared with the previous one
105 k(:,7) = feval (f, t + dt, varargout{2}, options.vfunarguments{:}); 110 k(:,7) = feval (f, t + dt, varargout{2}, args{:});
106 ## x_est according to Shampine 1986: 111 cc_prime = dt * c_prime;
107 ## varargout{3} = x + dt * ((1951/21600)*k(:,1) + (22642/50085)*k(:,3) 112 varargout{3} = x + k * cc_prime(:); # x_est
108 ## + (451/720)*k(:,4) - (12231/42400)*k(:,5)
109 ## + (649/6300)*k(:,6) + (1/60)*k(:,7));
110 varargout{3} = x + dt * ((5179/57600)*k(:,1) + (7571/16695)*k(:,3)
111 + (393/640)*k(:,4) - (92097/339200)*k(:,5)
112 + (187/2100)*k(:,6) + (1/40)*k(:,7)); # x_est
113 varargout{4} = k; 113 varargout{4} = k;
114 endif 114 endif
115 115
116 endfunction 116 endfunction
117 117