Mercurial > octave-nkf

comparison scripts/ode/private/integrate_adaptive.m @ 20637:756b052037fb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
avoid stepping beyond end of thspan in ode solvers * scripts/ode/private/integrate_adaptive.m: make sure not to step beyound end of tspan.
author Carlo de Falco <carlo.defalco@polimi.it>
date Sun, 11 Oct 2015 20:39:11 +0200
parents 43822bda4f65
children a260a6acb70f
comparison
equal deleted inserted replaced
20636:43822bda4f65 20637:756b052037fb
64 64
65 fixed_times = numel (tspan) > 2; 65 fixed_times = numel (tspan) > 2;
66 66
67 t_new = t_old = t = tspan(1); 67 t_new = t_old = t = tspan(1);
68 x_new = x_old = x = x0(:); 68 x_new = x_old = x = x0(:);
69 69
70 ## get first initial timestep 70 ## get first initial timestep
71 dt = starting_stepsize (order, func, t, x, options.AbsTol, 71 dt = starting_stepsize (order, func, t, x, options.AbsTol,
72 options.RelTol, options.vnormcontrol); 72 options.RelTol, options.vnormcontrol);
73 dt = odeget (options, "InitialStep", dt, "fast_not_empty"); 73 dt = odeget (options, "InitialStep", dt, "fast_not_empty");
74 74
75 dir = odeget (options, "vdirection", [], "fast"); 75 dir = odeget (options, "vdirection", [], "fast");
76 dt = dir * min (abs (dt), options.MaxStep); 76 dt = dir * min (abs (dt), options.MaxStep);
77 77
78 options.comp = 0.0; 78 options.comp = 0.0;
79 79
80 ## Factor multiplying the stepsize guess 80 ## Factor multiplying the stepsize guess
81 facmin = 0.8; 81 facmin = 0.8;
82 facmax = 1.5; 82 facmax = 1.5;
83 fac = 0.38^(1/(order+1)); # formula taken from Hairer 83 fac = 0.38^(1/(order+1)); # formula taken from Hairer
84 84
85 ## Initialize the OutputFcn 85 ## Initialize the OutputFcn
86 if (options.vhaveoutputfunction) 86 if (options.vhaveoutputfunction)
87 if (options.vhaveoutputselection) 87 if (options.vhaveoutputselection)
88 solution.vretout = x(options.OutputSel,end); 88 solution.vretout = x(options.OutputSel,end);
89 else 89 else
90 solution.vretout = x; 90 solution.vretout = x;
91 endif 91 endif
92 feval (options.OutputFcn, tspan, solution.vretout, 92 feval (options.OutputFcn, tspan, solution.vretout,
93 "init", options.vfunarguments{:}); 93 "init", options.vfunarguments{:});
94 endif 94 endif
100 endif 100 endif
101 101
102 if (options.vhavenonnegative) 102 if (options.vhavenonnegative)
103 nn = options.NonNegative; 103 nn = options.NonNegative;
104 endif 104 endif
105 105
106 solution.vcntloop = 2; 106 solution.vcntloop = 2;
107 solution.vcntcycles = 1; 107 solution.vcntcycles = 1;
108 solution.vcntsave = 2; 108 solution.vcntsave = 2;
109 solution.vunhandledtermination = true; 109 solution.vunhandledtermination = true;
110 ireject = 0; 110 ireject = 0;
111 111
112 k_vals = []; 112 k_vals = [];
113 iout = istep = 1; 113 iout = istep = 1;
114 while (dir * t_old < dir * tspan(end)) 114 while (dir * t_old < dir * tspan(end))
115 115
116 ## Compute integration step from t_old to t_new = t_old + dt 116 ## Compute integration step from t_old to t_new = t_old + dt
117 [t_new, options.comp] = kahan (t_old, options.comp, dt); 117 [t_new, options.comp] = kahan (t_old, options.comp, dt);
126 y_est(nn,end) = abs (y_est(nn,end)); 126 y_est(nn,end) = abs (y_est(nn,end));
127 endif 127 endif
128 128
129 err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol, 129 err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol,
130 options.vnormcontrol, x_est); 130 options.vnormcontrol, x_est);
131 131
132 ## Accepted solution only if err <= 1.0 132 ## Accepted solution only if err <= 1.0
133 if (err <= 1) 133 if (err <= 1)
134 134
135 solution.vcntloop++; 135 solution.vcntloop++;
136 ireject = 0; 136 ireject = 0;
137 137
138 ## if output time steps are fixed 138 ## if output time steps are fixed
139 if (fixed_times) 139 if (fixed_times)
140 140
141 t_caught = find ((tspan(iout:end) > t_old) 141 t_caught = find ((tspan(iout:end) > t_old)
142 & (tspan(iout:end) <= t_new)); 142 & (tspan(iout:end) <= t_new));
143 t_caught = t_caught + iout - 1; 143 t_caught = t_caught + iout - 1;
144 144
145 if (! isempty (t_caught)) 145 if (! isempty (t_caught))
146 t(t_caught) = tspan(t_caught); 146 t(t_caught) = tspan(t_caught);
147 iout = max (t_caught); 147 iout = max (t_caught);
148 x(:, t_caught) = ... 148 x(:, t_caught) = ...
149 ode_rk_interpolate (order, [t_old t_new], [x_old x_new], 149 ode_rk_interpolate (order, [t_old t_new], [x_old x_new],
166 && solution.vevent{1} == 1) 166 && solution.vevent{1} == 1)
167 t(id) = solution.vevent{3}(end); 167 t(id) = solution.vevent{3}(end);
168 t = t(1:id); 168 t = t(1:id);
169 x(:, id) = solution.vevent{4}(end, :).'; 169 x(:, id) = solution.vevent{4}(end, :).';
170 x = x(:,1:id); 170 x = x(:,1:id);
171 solution.vunhandledtermination = false; 171 solution.vunhandledtermination = false;
172 break_loop = true; 172 break_loop = true;
173 break; 173 break;
174 endif 174 endif
175 endfor 175 endfor
176 if (break_loop) 176 if (break_loop)
196 break; 196 break;
197 endif 197 endif
198 endif 198 endif
199 199
200 endif 200 endif
201 201
202 else 202 else
203 203
204 t(++istep) = t_new; 204 t(++istep) = t_new;
205 x(:, istep) = x_new; 205 x(:, istep) = x_new;
206 iout = istep; 206 iout = istep;
213 options.vfunarguments{:}); 213 options.vfunarguments{:});
214 if (! isempty (solution.vevent{1}) 214 if (! isempty (solution.vevent{1})
215 && solution.vevent{1} == 1) 215 && solution.vevent{1} == 1)
216 t(istep) = solution.vevent{3}(end); 216 t(istep) = solution.vevent{3}(end);
217 x(:, istep) = solution.vevent{4}(end, :).'; 217 x(:, istep) = solution.vevent{4}(end, :).';
218 solution.vunhandledtermination = false; 218 solution.vunhandledtermination = false;
219 break; 219 break;
220 endif 220 endif
221 endif 221 endif
222 222
223 ## Call plot. Stop integration if plot function 223 ## Call plot. Stop integration if plot function
246 x_old = x_new; 246 x_old = x_new;
247 k_vals = new_k_vals; 247 k_vals = new_k_vals;
248 248
249 solution.vcntloop = solution.vcntloop + 1; 249 solution.vcntloop = solution.vcntloop + 1;
250 vcntiter = 0; 250 vcntiter = 0;
251 251
252 else 252 else
253 253
254 ireject++; 254 ireject++;
255 255
256 ## Stop solving because in the last 5000 steps no successful valid 256 ## Stop solving because in the last 5000 steps no successful valid
265 " 'MaxStep' with the command 'odeset'.\n"], 265 " 'MaxStep' with the command 'odeset'.\n"],
266 t_old, tspan(end)); 266 t_old, tspan(end));
267 endif 267 endif
268 268
269 endif 269 endif
270 270
271 ## Compute next timestep, formula taken from Hairer 271 ## Compute next timestep, formula taken from Hairer
272 err += eps; # avoid divisions by zero 272 err += eps; # avoid divisions by zero
273 dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1)))); 273 dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1))));
274 dt = dir * min (abs (dt), options.MaxStep); 274 dt = dir * min (abs (dt), options.MaxStep);
275
276 ## make sure we don't go past tpan(end)
277 dt = dir * min (abs (dt), abs (tspan(end) - t_old));
275 278
276 endwhile 279 endwhile
277 280
281
278 ## Check if integration of the ode has been successful 282 ## Check if integration of the ode has been successful
279 if (dir * t(end) < dir * tspan(end)) 283 if (dir * t(end) < dir * tspan(end))
280 if (solution.vunhandledtermination == true) 284 if (solution.vunhandledtermination == true)
281 error ("integrate_adaptive:unexpected_termination", 285 error ("integrate_adaptive:unexpected_termination",
282 [" Solving has not been successful. The iterative", 286 [" Solving has not been successful. The iterative",
298 endif 302 endif
299 303
300 ## Set up return structure 304 ## Set up return structure
301 solution.t = t(:); 305 solution.t = t(:);
302 solution.x = x.'; 306 solution.x = x.';
303 307
304 endfunction 308 endfunction
305 309