Mercurial > hg > octave-nkf
comparison src/DLD-FUNCTIONS/dasrt.cc @ 4115:fc2048d4cd21
[project @ 2002-10-22 21:28:42 by jwe]
author | jwe |
---|---|
date | Tue, 22 Oct 2002 21:28:42 +0000 |
parents | c2562b2593e2 |
children | 944b276d8856 |
comparison
equal
deleted
inserted
replaced
4114:a32457362437 | 4115:fc2048d4cd21 |
---|---|
41 #include "utils.h" | 41 #include "utils.h" |
42 #include "variables.h" | 42 #include "variables.h" |
43 | 43 |
44 #include "DASRT-opts.cc" | 44 #include "DASRT-opts.cc" |
45 | 45 |
46 // Global pointers for user defined function required by dassl. | 46 // Global pointers for user defined function required by dasrt. |
47 static octave_function *dasrt_f; | 47 static octave_function *dasrt_f; |
48 static octave_function *dasrt_j; | 48 static octave_function *dasrt_j; |
49 static octave_function *dasrt_cf; | 49 static octave_function *dasrt_cf; |
50 | 50 |
51 // Is this a recursive call? | 51 // Is this a recursive call? |
233 } \ | 233 } \ |
234 while (0) | 234 while (0) |
235 | 235 |
236 DEFUN_DLD (dasrt, args, nargout, | 236 DEFUN_DLD (dasrt, args, nargout, |
237 "-*- texinfo -*-\n\ | 237 "-*- texinfo -*-\n\ |
238 @deftypefn {Loadable Function} {[@var{x}, @var{xdot}, @var{t}] =} dasrt (@var{fj} [, @var{g}], @var{x_0}, @var{xdot_0}, @var{t_out} [, @var{t_crit}])\n\ | 238 @deftypefn {Loadable Function} {[@var{x}, @var{xdot}, @var{t_out}, @var{istat}, @var{msg}] =} dasrt (@var{fcn} [, @var{g}], @var{x_0}, @var{xdot_0}, @var{t} [, @var{t_crit}])\n\ |
239 Solve a system of differential/algebraic equations with functional\n\ | 239 Solve the set of differential-algebraic equations\n\ |
240 stopping criteria.\n\ | 240 @tex\n\ |
241 \n\ | 241 $$ 0 = f (\\dot{x}, x, t) $$\n\ |
242 The function to be integrated must be of the form:\n\ | 242 with\n\ |
243 $$ x(t_0) = x_0, \\dot{x}(t_0) = \\dot{x}_0 $$\n\ | |
244 @end tex\n\ | |
245 @ifinfo\n\ | |
246 \n\ | |
243 @example\n\ | 247 @example\n\ |
244 @var{res} = f (@var{x}, @var{xdot}, @var{t}) = 0\n\ | 248 0 = f (xdot, x, t)\n\ |
245 @end example\n\ | 249 @end example\n\ |
246 \n\ | 250 \n\ |
247 The stopping condition must be of the form:\n\ | 251 with\n\ |
248 \n\ | 252 \n\ |
249 @example\n\ | 253 @example\n\ |
250 @var{res} = g (@var{x}, @var{t}) = 0\n\ | 254 x(t_0) = x_0, xdot(t_0) = xdot_0\n\ |
251 @end example\n\ | 255 @end example\n\ |
252 \n\ | 256 \n\ |
253 The Jacobian (if included) must be of the form:\n\ | 257 @end ifinfo\n\ |
258 with functional stopping criteria (root solving).\n\ | |
259 \n\ | |
260 The solution is returned in the matrices @var{x} and @var{xdot},\n\ | |
261 with each row in the result matrices corresponding to one of the\n\ | |
262 elements in the vector @var{t_out}. The first element of @var{t}\n\ | |
263 should be @math{t_0} and correspond to the initial state of the\n\ | |
264 system @var{x_0} and its derivative @var{xdot_0}, so that the first\n\ | |
265 row of the output @var{x} is @var{x_0} and the first row\n\ | |
266 of the output @var{xdot} is @var{xdot_0}.\n\ | |
267 \n\ | |
268 The vector @var{t} provides an upper limit on the length of the\n\ | |
269 integration. If the stopping condition is met, the vector\n\ | |
270 @var{t_out} will be shorter than @var{t}, and the final element of\n\ | |
271 @var{t_out} will be the point at which the stopping condition was met,\n\ | |
272 and may not correspond to any element of the vector @var{t}.\n\ | |
273 \n\ | |
274 The first argument, @var{fcn}, is a string that names the function to\n\ | |
275 call to compute the vector of residuals for the set of equations.\n\ | |
276 It must have the form\n\ | |
254 \n\ | 277 \n\ |
255 @example\n\ | 278 @example\n\ |
256 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{cj})\n\ | 279 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\ |
257 = df/dx + cj*df/dxdot\n\ | |
258 @end example\n\ | 280 @end example\n\ |
259 \n\ | 281 \n\ |
260 @noindent\n\ | 282 @noindent\n\ |
261 The following inputs are entered:\n\ | 283 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\ |
262 \n\ | 284 scalar.\n\ |
263 @table @var\n\ | 285 \n\ |
264 @item fj\n\ | 286 If @var{fcn} is a two-element string array, the first element names\n\ |
265 The string vector containing either @var{f} or both @var{f} and @var{j}.\n\ | 287 the function @math{f} described above, and the second element names\n\ |
266 \n\ | 288 a function to compute the modified Jacobian |
267 @item f\n\ | 289 \n\ |
268 The function to be integrated.\n\ | 290 @tex\n\ |
269 \n\ | 291 $$\n\ |
270 @item g\n\ | 292 J = {\\partial f \\over \\partial x}\n\ |
271 The function with the stopping conditions.\n\ | 293 + c {\\partial f \\over \\partial \\dot{x}}\n\ |
272 \n\ | 294 $$\n\ |
273 @item j\n\ | 295 @end tex\n\ |
274 The optional Jacobian function. If not included, it will be approximated\n\ | 296 @ifinfo\n\ |
275 by finite differences.\n\ | 297 \n\ |
276 \n\ | 298 @example\n\ |
277 @item x_0\n\ | 299 df df\n\ |
278 The initial state.\n\ | 300 jac = -- + c ------\n\ |
279 \n\ | 301 dx d xdot\n\ |
280 @item xdot_0\n\ | 302 @end example\n\ |
281 The time derivative of the initial state.\n\ | 303 \n\ |
282 \n\ | 304 @end ifinfo\n\ |
283 @item t_out\n\ | 305 \n\ |
284 The times at which the solution should be returned. This vector should\n\ | 306 The modified Jacobian function must have the form\n\ |
285 include the initial time a the first value.\n\ | 307 \n\ |
286 \n\ | 308 @example\n\ |
287 @item t_crit\n\ | 309 \n\ |
288 The times at which the function is non-smooth or poorly behaved.\n\ | 310 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\ |
289 \n\ | 311 \n\ |
290 @end table\n\ | 312 @end example\n\ |
291 \n\ | 313 \n\ |
292 @noindent\n\ | 314 The optional second argument names a function that defines the\n\ |
293 The following outputs are returned:\n\ | 315 constraint functions whose roots are desired during the integration.\n\ |
294 \n\ | 316 This function must have the form\n\ |
295 @table @var\n\ | 317 \n\ |
296 @item x\n\ | 318 @example\n\ |
297 The states at each time in @var{t}.\n\ | 319 @var{g_out} = g (@var{x}, @var{t})\n\ |
298 \n\ | 320 @end example\n\ |
299 @item xdot\n\ | 321 \n\ |
300 The time derivative of the states at each time in @var{t}.\n\ | 322 and return a vector of the constraint function values.\n\ |
301 \n\ | 323 If the value of any of the constraint functions changes sign, @sc{Dasrt}\n\ |
302 @item t\n\ | 324 will attempt to stop the integration at the point of the sign change.\n\ |
303 All the times in the requested output time vector up to the stopping\n\ | 325 \n\ |
304 criteria. The time at which the stopping criteria is achieved is returned\n\ | 326 If the name of the constraint function is omitted, @code{dasrt} solves\n\ |
305 as the last time in the vector.\n\ | 327 the saem problem as @code{daspk} or @code{dassl}. |
306 @end table\n\ | 328 \n\ |
329 Note that because of numerical errors in the constraint functions\n\ | |
330 due to roundoff and integration error, @sc{Dasrt} may return false\n\ | |
331 roots, or return the same root at two or more nearly equal values of\n\ | |
332 @var{T}. If such false roots are suspected, the user should consider\n\ | |
333 smaller error tolerances or higher precision in the evaluation of the\n\ | |
334 constraint functions.\n\ | |
335 \n\ | |
336 If a root of some constraint function defines the end of the problem,\n\ | |
337 the input to @sc{Dasrt} should nevertheless allow integration to a\n\ | |
338 point slightly past that root, so that @sc{Dasrt} can locate the root\n\ | |
339 by interpolation.\n\ | |
340 \n\ | |
341 The third and fourth arguments to @code{dasrt} specify the initial\n\ | |
342 condition of the states and their derivatives, and the fourth argument\n\ | |
343 specifies a vector of output times at which the solution is desired,\n\ | |
344 including the time corresponding to the initial condition.\n\ | |
345 \n\ | |
346 The set of initial states and derivatives are not strictly required to\n\ | |
347 be consistent. In practice, however, @sc{Dassl} is not very good at\n\ | |
348 determining a consistent set for you, so it is best if you ensure that\n\ | |
349 the initial values result in the function evaluating to zero.\n\ | |
350 \n\ | |
351 The sixth argument is optional, and may be used to specify a set of\n\ | |
352 times that the DAE solver should not integrate past. It is useful for\n\ | |
353 avoiding difficulties with singularities and points where there is a\n\ | |
354 discontinuity in the derivative.\n\ | |
355 \n\ | |
356 After a successful computation, the value of @var{istate} will be\n\ | |
357 greater than zero (consistent with the Fortran version of @sc{Dassl}).\n\ | |
358 \n\ | |
359 If the computation is not successful, the value of @var{istate} will be\n\ | |
360 less than zero and @var{msg} will contain additional information.\n\ | |
307 \n\ | 361 \n\ |
308 You can use the function @code{dasrt_options} to set optional\n\ | 362 You can use the function @code{dasrt_options} to set optional\n\ |
309 parameters for @code{dasrt}.\n\ | 363 parameters for @code{dasrt}.\n\ |
310 @end deftypefn") | 364 @end deftypefn\n\ |
365 @seealso{daspk, dasrt, lsode, odessa}") | |
311 { | 366 { |
312 octave_value_list retval; | 367 octave_value_list retval; |
313 | 368 |
314 unwind_protect::begin_frame ("Fdasrt"); | 369 unwind_protect::begin_frame ("Fdasrt"); |
315 | 370 |