Mercurial > octave-nkf

comparison scripts/ode/private/integrate_n_steps.m @ 20584:eb9e2d187ed2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use Octave coding conventions in scripts/ode/private dir. * AbsRel_Norm.m, fuzzy_compare.m, hermite_quartic_interpolation.m, integrate_adaptive.m, integrate_const.m, integrate_n_steps.m, kahan.m, ode_struct_value_check.m, odepkg_event_handle.m, odepkg_structure_check.m, runge_kutta_45_dorpri.m, starting_stepsize.m: Wrap long lines to < 80 chars. Use double quotes rather than single quotes where possible. Use ';' at end of keywords "return;" and "break;" Use '##" for stand-alone comments and '#' for end-of-line comments. Use two spaces after period before starting new sentence. Use '!' instead of '~' for logical negation. Use specific form of end (endif, endfor, etc.). Don't use line continuation marker '...' unless necessary.
author Rik <rik@octave.org>
date Sun, 04 Oct 2015 22:18:54 -0700
parents 25623ef2ff4f
children b7ac1e94266e
comparison
equal deleted inserted replaced
20583:d746695bf494 20584:eb9e2d187ed2
62 ## 62 ##
63 ## @seealso{integrate_adaptive, integrate_const} 63 ## @seealso{integrate_adaptive, integrate_const}
64 64
65 function solution = integrate_n_steps (stepper, func, t0, x0, dt, n, options) 65 function solution = integrate_n_steps (stepper, func, t0, x0, dt, n, options)
66 66
67 solution = struct; 67 solution = struct ();
68 68
69 ## first values for time and solution 69 ## first values for time and solution
70 x = x0(:); 70 x = x0(:);
71 t = t0; 71 t = t0;
72 72
73 vdirection = odeget (options, "vdirection", [], "fast"); 73 vdirection = odeget (options, "vdirection", [], "fast");
74 if (sign (dt) != vdirection) 74 if (sign (dt) != vdirection)
75 error("OdePkg:InvalidArgument", 75 error ("OdePkg:InvalidArgument", "option 'InitialStep' has a wrong sign");
76 "option ''InitialStep'' has a wrong sign");
77 endif 76 endif
78 77
79 comp = 0.0; 78 comp = 0.0;
80 tk = t0; 79 tk = t0;
81 options.comp = comp; 80 options.comp = comp;
137 endif 136 endif
138 solution.vcntloop = solution.vcntloop + 1; 137 solution.vcntloop = solution.vcntloop + 1;
139 vcntiter = 0; 138 vcntiter = 0;
140 139
141 ## Call plot only if a valid result has been found, therefore this code 140 ## Call plot only if a valid result has been found, therefore this code
142 ## fragment has moved here. Stop integration if plot function returns false 141 ## fragment has moved here. Stop integration if plot function returns false
143 if (options.vhaveoutputfunction) 142 if (options.vhaveoutputfunction)
144 for vcnt = 0:options.Refine # Approximation between told and t 143 for vcnt = 0:options.Refine # Approximation between told and t
145 if (options.vhaverefine) # Do interpolation 144 if (options.vhaverefine) # Do interpolation
146 vapproxtime = (vcnt + 1) / (options.Refine + 2); 145 vapproxtime = (vcnt + 1) / (options.Refine + 2);
147 vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ... 146 vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ...
155 vapproxvals = vapproxvals(options.OutputSel); 154 vapproxvals = vapproxvals(options.OutputSel);
156 endif 155 endif
157 vpltret = feval (options.OutputFcn, vapproxtime, vapproxvals, [], 156 vpltret = feval (options.OutputFcn, vapproxtime, vapproxvals, [],
158 options.vfunarguments{:}); 157 options.vfunarguments{:});
159 if (vpltret) # Leave refinement loop 158 if (vpltret) # Leave refinement loop
160 break 159 break;
161 endif 160 endif
162 endfor 161 endfor
163 if (vpltret) # Leave main loop 162 if (vpltret) # Leave main loop
164 solution.vunhandledtermination = false; 163 solution.vunhandledtermination = false;
165 break 164 break;
166 endif 165 endif
167 endif 166 endif
168 167
169 ## Call event only if a valid result has been found, therefore this 168 ## Call event only if a valid result has been found, therefore this
170 ## code fragment has moved here. Stop integration if veventbreak is 169 ## code fragment has moved here. Stop integration if veventbreak is
171 ## true 170 ## true
172 if (options.vhaveeventfunction) 171 if (options.vhaveeventfunction)
173 solution.vevent = odepkg_event_handle (options.Events, t(end), x(:,end), 172 solution.vevent = odepkg_event_handle (options.Events, t(end), x(:,end),
174 [], options.vfunarguments{:}); 173 [], options.vfunarguments{:});
175 if (! isempty (solution.vevent{1}) 174 if (! isempty (solution.vevent{1})
176 && solution.vevent{1} == 1) 175 && solution.vevent{1} == 1)
177 t(solution.vcntloop-1,:) = solution.vevent{3}(end,:); 176 t(solution.vcntloop-1,:) = solution.vevent{3}(end,:);
178 x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)'; 177 x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)';
179 solution.vunhandledtermination = false; 178 solution.vunhandledtermination = false;
180 break 179 break;
181 endif 180 endif
182 endif 181 endif
183 182
184 ## Update counters that count the number of iteration cycles 183 ## Update counters that count the number of iteration cycles
185 solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics 184 solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics
186 vcntiter = vcntiter + 1; # Needed to find iteration problems 185 vcntiter = vcntiter + 1; # Needed to find iteration problems
187 186
188 ## Stop solving because the last 1000 steps no successful valid 187 ## Stop solving because the last 1000 steps no successful valid
189 ## value has been found 188 ## value has been found
190 if (vcntiter >= 5000) 189 if (vcntiter >= 5000)
191 error (["Solving has not been successful. The iterative", 190 error (["Solving has not been successful. The iterative",
192 " integration loop exited at time t = %f before endpoint at", 191 " integration loop exited at time t = %f before endpoint at",
193 " tend = %f was reached. This happened because the iterative", 192 " tend = %f was reached. This happened because the iterative",
194 " integration loop does not find a valid solution at this time", 193 " integration loop does not find a valid solution at this time",
195 " stamp. Try to reduce the value of ''InitialStep'' and/or", 194 " stamp. Try to reduce the value of 'InitialStep' and/or",
196 " ''MaxStep'' with the command ''odeset''.\n"], 195 " 'MaxStep' with the command 'odeset'.\n"],
197 s(end), tspan(end)); 196 s(end), tspan(end));
198 endif 197 endif
199 endfor 198 endfor
200 199
201 ## Check if integration of the ode has been successful 200 ## Check if integration of the ode has been successful
202 #if (vdirection * z(end) < vdirection * tspan(end)) 201 #if (vdirection * z(end) < vdirection * tspan(end))
203 # if (solution.vunhandledtermination == true) 202 # if (solution.vunhandledtermination == true)
204 # error ("OdePkg:InvalidArgument", 203 # error ("OdePkg:InvalidArgument",
205 # ["Solving has not been successful. The iterative", 204 # ["Solving has not been successful. The iterative",
206 # " integration loop exited at time t = %f", 205 # " integration loop exited at time t = %f",
207 # " before endpoint at tend = %f was reached. This may", 206 # " before endpoint at tend = %f was reached. This may",
208 # " happen if the stepsize grows smaller than defined in", 207 # " happen if the stepsize grows smaller than defined in",
209 # " vminstepsize. Try to reduce the value of ''InitialStep''", 208 # " vminstepsize. Try to reduce the value of 'InitialStep'",
210 # " and/or ''MaxStep'' with the command ''odeset''.\n"], 209 # " and/or 'MaxStep' with the command 'odeset'.\n"],
211 # z(end), tspan(end)); 210 # z(end), tspan(end));
212 # else 211 # else
213 # warning ("OdePkg:InvalidArgument", 212 # warning ("OdePkg:InvalidArgument",
214 # ["Solver has been stopped by a call of ''break'' in the main", 213 # ["Solver has been stopped by a call of 'break' in the main",
215 # " iteration loop at time t = %f before endpoint at tend = %f", 214 # " iteration loop at time t = %f before endpoint at tend = %f",
216 # " was reached. This may happen because the @odeplot function", 215 # " was reached. This may happen because the @odeplot function",
217 # " returned ''true'' or the @event function returned", 216 # " returned 'true' or the @event function returned",
218 # " ''true''.\n"], 217 # " 'true'.\n"],
219 # z(end), tspan(end)); 218 # z(end), tspan(end));
220 # endif 219 # endif
221 #endif 220 #endif
222 221
223 solution.t = t; 222 solution.t = t;
224 solution.x = x'; 223 solution.x = x';
225 224
226 endfunction 225 endfunction
227 226
228 ## Local Variables: ***
229 ## mode: octave ***
230 ## End: ***