Mercurial > octave-nkf

comparison scripts/ode/private/integrate_adaptive.m @ 20568:fcb792acab9b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Moving ode45, odeset, odeget, and levenshtein from odepkg to core. * libinterp/corefcn/levenshtein.cc: move function from odepkg into core * libinterp/corefcn/module.mk: include levenshtein.cc * scripts/ode: move ode45, odeset, odeget, and all dependencies from odepkg into core * scripts/module.mk: include them * doc/interpreter/diffeq.txi: add documentation for ode45, odeset, odeget * NEWS: announce functions included with this changeset * scripts/help/__unimplemented__.m: removed new functions
author jcorno <jacopo.corno@gmail.com>
date Thu, 24 Sep 2015 12:58:46 +0200
parents
children 6256f6e366ac
comparison
equal deleted inserted replaced
20567:2480bbcd1333 20568:fcb792acab9b
1 ## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it>
2 ## OdePkg - A package for solving ordinary differential equations and more
3 ##
4 ## This program is free software; you can redistribute it and/or modify
5 ## it under the terms of the GNU General Public License as published by
6 ## the Free Software Foundation; either version 2 of the License, or
7 ## (at your option) any later version.
8 ##
9 ## This program is distributed in the hope that it will be useful,
10 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ## GNU General Public License for more details.
13 ##
14 ## You should have received a copy of the GNU General Public License
15 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
16
17
18 ## -*- texinfo -*-
19 ## @deftypefn {Command} {[@var{t}, @var{y}] =}
20 ## integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@fun}, @var{tspan},
21 ## @var{x0}, @var{options})
22 ##
23 ## This function file can be called by a ODE solver function in order to
24 ## integrate the set of ODEs on the interval @var{[t0,t1]} with an
25 ## adaptive timestep.
26 ##
27 ## This function must be called with two output arguments: @var{t} and @var{y}.
28 ## Variable @var{t} is a column vector and contains the time stamps, instead
29 ## @var{y} is a matrix in which each column refers to a different unknown
30 ## of the problem and the rows number is the same of @var{t} rows number so
31 ## that each row of @var{y} contains the values of all unknowns at the time
32 ## value contained in the corresponding row in @var{t}.
33 ##
34 ## The first input argument must be a function_handle or an inline function
35 ## representing the stepper, that is the function responsible for step-by-step
36 ## integration. This function discriminates one method from the others.
37 ##
38 ## The second input argument is the order of the stepper. It is needed
39 ## to compute the adaptive timesteps.
40 ##
41 ## The third input argument is a function_handle or an inline function that
42 ## defines the set of ODE:
43 ## @ifhtml
44 ## @example
45 ## @math{y' = f(t,y)}
46 ## @end example
47 ## @end ifhtml
48 ## @ifnothtml
49 ## @math{y' = f(t,y)}.
50 ## @end ifnothtml
51 ##
52 ## The fourth input argument is the time vector which defines integration
53 ## interval, that is @var{[tspan(1),tspan(end)]} and all the intermediate
54 ## elements are taken as times at which the solution is required.
55 ##
56 ## The fifth argument represents the initial conditions for the ODEs and the
57 ## last input argument contains some options that may be needed for the stepper.
58 ##
59 ## @end deftypefn
60 ##
61 ## @seealso{integrate_const, integrate_n_steps}
62
63 function solution = integrate_adaptive (stepper, order, func, tspan, x0, options)
64
65 solution = struct;
66
67 ## first values for time and solution
68 t = tspan(1);
69 x = x0(:);
70
71 ## get first initial timestep
72 dt = odeget (options, "InitialStep",
73 starting_stepsize (order, func, t, x, options.AbsTol,
74 options.RelTol, options.vnormcontrol),
75 "fast_not_empty");
76 vdirection = odeget (options, "vdirection", [], "fast");
77 if (sign (dt) != vdirection)
78 dt = -dt;
79 endif
80 dt = vdirection * min (abs (dt), options.MaxStep);
81
82 ## set parameters
83 k = length (tspan);
84 counter = 2;
85 comp = 0.0;
86 tk = tspan(1);
87 options.comp = comp;
88
89 ## factor multiplying the stepsize guess
90 facmin = 0.8;
91 fac = 0.38^(1/(order+1)); ## formula taken from Hairer
92 t_caught = false;
93
94
95 ## Initialize the OutputFcn
96 if (options.vhaveoutputfunction)
97 if (options.vhaveoutputselection)
98 solution.vretout = x(options.OutputSel,end);
99 else
100 solution.vretout = x;
101 endif
102 feval (options.OutputFcn, tspan, solution.vretout,
103 "init", options.vfunarguments{:});
104 endif
105
106 ## Initialize the EventFcn
107 if (options.vhaveeventfunction)
108 odepkg_event_handle (options.Events, t(end), x,
109 "init", options.vfunarguments{:});
110 endif
111
112 solution.vcntloop = 2;
113 solution.vcntcycles = 1;
114 vcntiter = 0;
115 solution.vunhandledtermination = true;
116 solution.vcntsave = 2;
117
118 z = t;
119 u = x;
120
121 k_vals = feval (func, t , x, options.vfunarguments{:});
122
123 while (counter <= k)
124 facmax = 1.5;
125
126 ## compute integration step from t to t+dt
127 [s, y, y_est, k_vals] = stepper (func, z(end), u(:,end),
128 dt, options, k_vals);
129
130 if (options.vhavenonnegative)
131 x(options.NonNegative,end) = abs (x(options.NonNegative,end));
132 y(options.NonNegative,end) = abs (y(options.NonNegative,end));
133 y_est(options.NonNegative,end) = abs (y_est(options.NonNegative,end));
134 endif
135
136 if (options.vhaveoutputfunction && options.vhaverefine)
137 vSaveVUForRefine = u(:,end);
138 endif
139
140 err = AbsRel_Norm (y(:,end), u(:,end), options.AbsTol, options.RelTol,
141 options.vnormcontrol, y_est(:,end));
142
143 ## solution accepted only if the error is less or equal to 1.0
144 if (err <= 1)
145
146 [tk, comp] = kahan (tk, comp, dt);
147 options.comp = comp;
148 s(end) = tk;
149
150 ## values on this interval for time and solution
151 z = [z(end);s];
152 u = [u(:,end),y];
153
154 ## if next tspan value is caught, update counter
155 if ((z(end) == tspan(counter))
156 || (abs (z(end) - tspan(counter)) /
157 (max (abs (z(end)), abs (tspan(counter)))) < 8*eps) )
158 counter++;
159
160 ## if there is an element in time vector at which the solution is required
161 ## the program must compute this solution before going on with next steps
162 elseif (vdirection * z(end) > vdirection * tspan(counter))
163 ## initialize counter for the following cycle
164 i = 2;
165 while (i <= length (z))
166
167 ## if next tspan value is caught, update counter
168 if ((counter <= k)
169 && ((z(i) == tspan(counter))
170 || (abs (z(i) - tspan(counter)) /
171 (max (abs (z(i)), abs (tspan(counter)))) < 8*eps)) )
172 counter++;
173 endif
174 ## else, loop until there are requested values inside this subinterval
175 while ((counter <= k)
176 && (vdirection * z(i) > vdirection * tspan(counter)))
177 ## choose interpolation scheme according to order of the solver
178 switch order
179 case 1
180 u_interp = linear_interpolation ([z(i-1) z(i)],
181 [u(:,i-1) u(:,i)],
182 tspan(counter));
183 case 2
184 if (! isempty (k_vals))
185 der = k_vals(:,1);
186 else
187 der = feval (func, z(i-1) , u(:,i-1),
188 options.vfunarguments{:});
189 endif
190 u_interp = quadratic_interpolation ([z(i-1) z(i)],
191 [u(:,i-1) u(:,i)],
192 der, tspan(counter));
193 case 3
194 u_interp = ...
195 hermite_cubic_interpolation ([z(i-1) z(i)],
196 [u(:,i-1) u(:,i)],
197 [k_vals(:,1) k_vals(:,end)],
198 tspan(counter));
199 case 4
200 ## if ode45 is used without local extrapolation this function
201 ## doesn't require a new function evaluation.
202 u_interp = dorpri_interpolation ([z(i-1) z(i)],
203 [u(:,i-1) u(:,i)],
204 k_vals, tspan(counter));
205 case 5
206 ## ode45 with Dormand-Prince scheme:
207 ## 4th order approximation of y in t+dt/2 as proposed by
208 ## Shampine in Lawrence, Shampine, "Some Practical
209 ## Runge-Kutta Formulas", 1986.
210 u_half = u(:,i-1) ...
211 + 1/2*dt*((6025192743/30085553152) * k_vals(:,1)
212 + (51252292925/65400821598) * k_vals(:,3)
213 - (2691868925/45128329728) * k_vals(:,4)
214 + (187940372067/1594534317056) * k_vals(:,5)
215 - (1776094331/19743644256) * k_vals(:,6)
216 + (11237099/235043384) * k_vals(:,7));
217 u_interp = ...
218 hermite_quartic_interpolation ([z(i-1) z(i)],
219 [u(:,i-1) u_half u(:,i)],
220 [k_vals(:,1) k_vals(:,end)],
221 tspan(counter));
222
223 ## it is also possible to do a new function evaluation and use
224 ## the quintic hermite interpolator
225 ## f_half = feval (func, t+1/2*dt, u_half,
226 ## options.vfunarguments{:});
227 ## u_interp =
228 ## hermite_quintic_interpolation ([z(i-1) z(i)],
229 ## [u(:,i-1) u_half u(:,i)],
230 ## [k_vals(:,1) f_half k_vals(:,end)],
231 ## tspan(counter));
232 otherwise
233 warning ("High order interpolation not yet implemented: ",
234 "using cubic iterpolation instead");
235 der(:,1) = feval (func, z(i-1) , u(:,i-1),
236 options.vfunarguments{:});
237 der(:,2) = feval (func, z(i) , u(:,i),
238 options.vfunarguments{:});
239 u_interp = ...
240 hermite_cubic_interpolation ([z(i-1) z(i)],
241 [u(:,i-1) u(:,i)],
242 der, tspan(counter));
243 endswitch
244
245 ## add the interpolated value of the solution
246 u = [u(:,1:i-1), u_interp, u(:,i:end)];
247
248 ## add the time requested
249 z = [z(1:i-1);tspan(counter);z(i:end)];
250
251 ## update counters
252 counter++;
253 i++;
254 endwhile
255
256 ## if new time requested is not out of this interval
257 if ((counter <= k)
258 && (vdirection * z(end) > vdirection * tspan(counter)))
259 ## update the counter
260 i++;
261 else
262 ## stop the cycle and go on with the next iteration
263 i = length (z) + 1;
264 endif
265
266 endwhile
267 endif
268
269 if (mod (solution.vcntloop-1, options.OutputSave) == 0)
270 x = [x,u(:,2:end)];
271 t = [t;z(2:end)];
272 solution.vcntsave = solution.vcntsave + 1;
273 endif
274 solution.vcntloop = solution.vcntloop + 1;
275 vcntiter = 0;
276
277 ## Call plot only if a valid result has been found, therefore this
278 ## code fragment has moved here. Stop integration if plot function
279 ## returns false
280 if (options.vhaveoutputfunction)
281 for vcnt = 0:options.Refine # Approximation between told and t
282 if (options.vhaverefine) # Do interpolation
283 vapproxtime = (vcnt + 1) / (options.Refine + 2);
284 vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ...
285 + (vapproxtime) * y(:,end);
286 vapproxtime = s(end) + vapproxtime * dt;
287 else
288 vapproxvals = x(:,end);
289 vapproxtime = t(end);
290 endif
291 if (options.vhaveoutputselection)
292 vapproxvals = vapproxvals(options.OutputSel);
293 endif
294 vpltret = feval (options.OutputFcn, vapproxtime,
295 vapproxvals, [], options.vfunarguments{:});
296 if (vpltret) # Leave refinement loop
297 break
298 endif
299 endfor
300 if (vpltret) # Leave main loop
301 solution.vunhandledtermination = false;
302 break
303 endif
304 endif
305
306 ## Call event only if a valid result has been found, therefore this
307 ## code fragment has moved here. Stop integration if veventbreak is
308 ## true
309 if (options.vhaveeventfunction)
310 solution.vevent = odepkg_event_handle (options.Events, t(end),
311 x(:,end), [], options.vfunarguments{:});
312 if (! isempty (solution.vevent{1})
313 && solution.vevent{1} == 1)
314 t(solution.vcntloop-1,:) = solution.vevent{3}(end,:);
315 x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)';
316 solution.vunhandledtermination = false;
317 break
318 endif
319 endif
320
321 else
322 facmax = 1.0;
323 endif
324
325 ## Compute next timestep, formula taken from Hairer
326 err += eps; # adding an eps to avoid divisions by zero
327 dt = dt * min (facmax, max (facmin,
328 fac * (1 / err)^(1 / (order + 1))));
329 dt = vdirection * min (abs (dt), options.MaxStep);
330
331 ## Update counters that count the number of iteration cycles
332 solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics
333 vcntiter = vcntiter + 1; # Needed to find iteration problems
334
335 ## Stop solving because in the last 1000 steps no successful valid
336 ## value has been found
337 if (vcntiter >= 5000)
338 error (["Solving has not been successful. The iterative",
339 " integration loop exited at time t = %f before endpoint at",
340 " tend = %f was reached. This happened because the iterative",
341 " integration loop does not find a valid solution at this time",
342 " stamp. Try to reduce the value of ''InitialStep'' and/or",
343 " ''MaxStep'' with the command ''odeset''.\n"],
344 s(end), tspan(end));
345 endif
346
347 ## if this is the last iteration, save the length of last interval
348 if (counter > k)
349 j = length (z);
350 endif
351 endwhile
352
353 ## Check if integration of the ode has been successful
354 if (vdirection * z(end) < vdirection * tspan(end))
355 if (solution.vunhandledtermination == true)
356 error ("OdePkg:InvalidArgument",
357 ["Solving has not been successful. The iterative",
358 " integration loop exited at time t = %f",
359 " before endpoint at tend = %f was reached. This may",
360 " happen if the stepsize grows smaller than defined in",
361 " vminstepsize. Try to reduce the value of ''InitialStep''",
362 " and/or ''MaxStep'' with the command ''odeset''.\n"],
363 z(end), tspan(end));
364 else
365 warning ("OdePkg:InvalidArgument",
366 ["Solver has been stopped by a call of ''break'' in the main",
367 " iteration loop at time t = %f before endpoint at tend = %f ",
368 " was reached. This may happen because the @odeplot function",
369 " returned ''true'' or the @event function returned",
370 " ''true''.\n"],
371 z(end), tspan(end));
372 endif
373 endif
374
375 ## Compute how many values are out of time inerval
376 d = vdirection * t((end-(j-1)):end) > vdirection * tspan(end)*ones (j, 1);
377 f = sum (d);
378
379 ## Remove not-requested values of time and solution
380 solution.t = t(1:end-f);
381 solution.x = x(:,1:end-f)';
382
383 endfunction
384
385 ## Local Variables: ***
386 ## mode: octave ***
387 ## End: ***