Mercurial > octave-nkf

comparison scripts/ode/private/integrate_adaptive.m @ 20639:a260a6acb70f

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix test failures introduced by a22d8a2eb0e5 * scripts/ode/private/integrate_adaptive.m: fix stepping backwards, fix invocation of OutputFcn, fix text of some error messages * scripts/ode/private/integrate_const.m: remove use of option OutputSave * scripts/ode/private/integrate_n_steps.m: remove use of option OutputSave
author Carlo de Falco <carlo.defalco@polimi.it>
date Sun, 11 Oct 2015 23:09:01 +0200
parents 756b052037fb
children
comparison
equal deleted inserted replaced
20638:d30fc2c11455 20639:a260a6acb70f
119 stepper (func, t_old, x_old, dt, options, k_vals, t_new); 119 stepper (func, t_old, x_old, dt, options, k_vals, t_new);
120 120
121 solution.vcntcycles++; 121 solution.vcntcycles++;
122 122
123 if (options.vhavenonnegative) 123 if (options.vhavenonnegative)
124 x(nn,end) = abs (x(nn,end)); 124 x_new(nn, end) = abs (x_new(nn, end));
125 y(nn,end) = abs (y(nn,end)); 125 x_est(nn, end) = abs (x_est(nn, end));
126 y_est(nn,end) = abs (y_est(nn,end));
127 endif 126 endif
128 127
129 err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol, 128 err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol,
130 options.vnormcontrol, x_est); 129 options.vnormcontrol, x_est);
131 130
136 ireject = 0; 135 ireject = 0;
137 136
138 ## if output time steps are fixed 137 ## if output time steps are fixed
139 if (fixed_times) 138 if (fixed_times)
140 139
141 t_caught = find ((tspan(iout:end) > t_old) 140 t_caught = find ((dir * tspan(iout:end) > dir * t_old)
142 & (tspan(iout:end) <= t_new)); 141 & (dir * tspan(iout:end) <= dir * t_new));
143 t_caught = t_caught + iout - 1; 142 t_caught = t_caught + iout - 1;
144 143
145 if (! isempty (t_caught)) 144 if (! isempty (t_caught))
146 t(t_caught) = tspan(t_caught); 145 t(t_caught) = tspan(t_caught);
147 iout = max (t_caught); 146 iout = max (t_caught);
182 ## returns true. 181 ## returns true.
183 if (options.vhaveoutputfunction) 182 if (options.vhaveoutputfunction)
184 vcnt = options.Refine + 1; 183 vcnt = options.Refine + 1;
185 vapproxtime = linspace (t_old, t_new, vcnt); 184 vapproxtime = linspace (t_old, t_new, vcnt);
186 vapproxvals = interp1 ([t_old, t(t_caught), t_new], 185 vapproxvals = interp1 ([t_old, t(t_caught), t_new],
187 [x_old, x(:, t_caught), x_new], 186 [x_old, x(:, t_caught), x_new] .',
188 vapproxtime, 'linear'); 187 vapproxtime, 'linear') .';
189 if (options.vhaveoutputselection) 188 if (options.vhaveoutputselection)
190 vapproxvals = vapproxvals(options.OutputSel); 189 vapproxvals = vapproxvals(options.OutputSel, :);
191 endif 190 endif
192 vpltret = feval (options.OutputFcn, vapproxtime, 191 for ii = 1:numel (vapproxtime)
193 vapproxvals, [], options.vfunarguments{:}); 192 vpltret = feval (options.OutputFcn, vapproxtime(ii),
193 vapproxvals(:, ii), [],
194 options.vfunarguments{:});
195 endfor
194 if (vpltret) # Leave main loop 196 if (vpltret) # Leave main loop
195 solution.vunhandledtermination = false; 197 solution.vunhandledtermination = false;
196 break; 198 break;
197 endif 199 endif
198 endif 200 endif
224 ## returns true. 226 ## returns true.
225 if (options.vhaveoutputfunction) 227 if (options.vhaveoutputfunction)
226 vcnt = options.Refine + 1; 228 vcnt = options.Refine + 1;
227 vapproxtime = linspace (t_old, t_new, vcnt); 229 vapproxtime = linspace (t_old, t_new, vcnt);
228 vapproxvals = interp1 ([t_old, t_new], 230 vapproxvals = interp1 ([t_old, t_new],
229 [x_old, x_new], 231 [x_old, x_new] .',
230 vapproxtime, 'linear'); 232 vapproxtime, 'linear') .';
231 if (options.vhaveoutputselection) 233 if (options.vhaveoutputselection)
232 vapproxvals = vapproxvals(options.OutputSel); 234 vapproxvals = vapproxvals(options.OutputSel, :);
233 endif 235 endif
234 vpltret = feval (options.OutputFcn, vapproxtime, 236 for ii = 1:numel (vapproxtime)
235 vapproxvals, [], options.vfunarguments{:}); 237 vpltret = feval (options.OutputFcn, vapproxtime(ii),
238 vapproxvals(:, ii), [], options.vfunarguments{:});
239 endfor
236 if (vpltret) # Leave main loop 240 if (vpltret) # Leave main loop
237 solution.vunhandledtermination = false; 241 solution.vunhandledtermination = false;
238 break; 242 break;
239 endif 243 endif
240 endif 244 endif
281 285
282 ## Check if integration of the ode has been successful 286 ## Check if integration of the ode has been successful
283 if (dir * t(end) < dir * tspan(end)) 287 if (dir * t(end) < dir * tspan(end))
284 if (solution.vunhandledtermination == true) 288 if (solution.vunhandledtermination == true)
285 error ("integrate_adaptive:unexpected_termination", 289 error ("integrate_adaptive:unexpected_termination",
286 [" Solving has not been successful. The iterative", 290 [" Solving has not been successful. The iterative", ...
287 " integration loop exited at time t = %f", 291 " integration loop exited at time t = %f ", ...
288 " before endpoint at tend = %f was reached. This may", 292 " before endpoint at tend = %f was reached. This may", ...
289 " happen if the stepsize grows too small. ", 293 " happen if the stepsize grows too small. ", ...
290 " Try to reduce the value of 'InitialStep'", 294 " Try to reduce the value of 'InitialStep'", ...
291 " and/or 'MaxStep' with the command 'odeset'.\n"], 295 " and/or 'MaxStep' with the command 'odeset'.\n"],
292 t(end), tspan(end)); 296 t(end), tspan(end));
293 else 297 else
294 warning ("integrate_adaptive:unexpected_termination", 298 warning ("integrate_adaptive:unexpected_termination",
295 ["Solver has been stopped by a call of 'break' in the main", 299 ["Solver has been stopped by a call of 'break' ", ...
296 " iteration loop at time t = %f before endpoint at tend = %f ", 300 " in the main iteration loop at time t = %f before", ...
297 " was reached. This may happen because the @odeplot function", 301 " endpoint at tend = %f was reached. This may", ...
298 " returned 'true' or the @event function returned", 302 " happen because the @odeplot function", ...
303 " returned 'true' or the @event function returned", ...
299 " 'true'.\n"], 304 " 'true'.\n"],
300 t(end), tspan(end)); 305 t(end), tspan(end));
301 endif 306 endif
302 endif 307 endif
303 308