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