GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
dasrt.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2002-2015 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include <iostream>
28 #include <string>
29 
30 #include "DASRT.h"
31 #include "lo-mappers.h"
32 
33 #include "defun.h"
34 #include "error.h"
35 #include "gripes.h"
36 #include "oct-obj.h"
37 #include "ov-fcn.h"
38 #include "ov-cell.h"
39 #include "pager.h"
40 #include "parse.h"
41 #include "unwind-prot.h"
42 #include "utils.h"
43 #include "variables.h"
44 
45 #include "DASRT-opts.cc"
46 
47 // Global pointers for user defined function required by dasrt.
51 
52 // Have we warned about imaginary values returned from user function?
53 static bool warned_fcn_imaginary = false;
54 static bool warned_jac_imaginary = false;
55 static bool warned_cf_imaginary = false;
56 
57 // Is this a recursive call?
58 static int call_depth = 0;
59 
60 static ColumnVector
61 dasrt_user_f (const ColumnVector& x, const ColumnVector& xdot,
62  double t, octave_idx_type&)
63 {
64  ColumnVector retval;
65 
66  assert (x.capacity () == xdot.capacity ());
67 
68  octave_value_list args;
69 
70  args(2) = t;
71  args(1) = xdot;
72  args(0) = x;
73 
74  if (dasrt_f)
75  {
76  octave_value_list tmp = dasrt_f->do_multi_index_op (1, args);
77 
78  if (error_state)
79  {
80  gripe_user_supplied_eval ("dasrt");
81  return retval;
82  }
83 
84  if (tmp.length () > 0 && tmp(0).is_defined ())
85  {
86  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
87  {
88  warning ("dasrt: ignoring imaginary part returned from user-supplied function");
89  warned_fcn_imaginary = true;
90  }
91 
92  retval = ColumnVector (tmp(0).vector_value ());
93 
94  if (error_state || retval.length () == 0)
95  gripe_user_supplied_eval ("dasrt");
96  }
97  else
98  gripe_user_supplied_eval ("dasrt");
99  }
100 
101  return retval;
102 }
103 
104 static ColumnVector
105 dasrt_user_cf (const ColumnVector& x, double t)
106 {
107  ColumnVector retval;
108 
109  octave_value_list args;
110 
111  args(1) = t;
112  args(0) = x;
113 
114  if (dasrt_cf)
115  {
116  octave_value_list tmp = dasrt_cf->do_multi_index_op (1, args);
117 
118  if (error_state)
119  {
120  gripe_user_supplied_eval ("dasrt");
121  return retval;
122  }
123 
124  if (tmp.length () > 0 && tmp(0).is_defined ())
125  {
126  if (! warned_cf_imaginary && tmp(0).is_complex_type ())
127  {
128  warning ("dasrt: ignoring imaginary part returned from user-supplied constraint function");
129  warned_cf_imaginary = true;
130  }
131 
132  retval = ColumnVector (tmp(0).vector_value ());
133 
134  if (error_state || retval.length () == 0)
135  gripe_user_supplied_eval ("dasrt");
136  }
137  else
138  gripe_user_supplied_eval ("dasrt");
139  }
140 
141  return retval;
142 }
143 
144 static Matrix
146  double t, double cj)
147 {
148  Matrix retval;
149 
150  assert (x.capacity () == xdot.capacity ());
151 
152  octave_value_list args;
153 
154  args(3) = cj;
155  args(2) = t;
156  args(1) = xdot;
157  args(0) = x;
158 
159  if (dasrt_j)
160  {
161  octave_value_list tmp = dasrt_j->do_multi_index_op (1, args);
162 
163  if (error_state)
164  {
165  gripe_user_supplied_eval ("dasrt");
166  return retval;
167  }
168 
169  int tlen = tmp.length ();
170  if (tlen > 0 && tmp(0).is_defined ())
171  {
172  if (! warned_jac_imaginary && tmp(0).is_complex_type ())
173  {
174  warning ("dasrt: ignoring imaginary part returned from user-supplied jacobian function");
175  warned_jac_imaginary = true;
176  }
177 
178  retval = tmp(0).matrix_value ();
179 
180  if (error_state || retval.length () == 0)
181  gripe_user_supplied_eval ("dasrt");
182  }
183  else
184  gripe_user_supplied_eval ("dasrt");
185  }
186 
187  return retval;
188 }
189 
190 #define DASRT_ABORT \
191  return retval
192 
193 #define DASRT_ABORT1(msg) \
194  do \
195  { \
196  ::error ("dasrt: " msg); \
197  DASRT_ABORT; \
198  } \
199  while (0)
200 
201 #define DASRT_ABORT2(fmt, arg) \
202  do \
203  { \
204  ::error ("dasrt: " fmt, arg); \
205  DASRT_ABORT; \
206  } \
207  while (0)
208 
209 DEFUN (dasrt, args, nargout,
210  "-*- texinfo -*-\n\
211 @deftypefn {Built-in Function} {[@var{x}, @var{xdot}, @var{t_out}, @var{istat}, @var{msg}] =} dasrt (@var{fcn}, [], @var{x_0}, @var{xdot_0}, @var{t})\n\
212 @deftypefnx {Built-in Function} {@dots{} =} dasrt (@var{fcn}, @var{g}, @var{x_0}, @var{xdot_0}, @var{t})\n\
213 @deftypefnx {Built-in Function} {@dots{} =} dasrt (@var{fcn}, [], @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})\n\
214 @deftypefnx {Built-in Function} {@dots{} =} dasrt (@var{fcn}, @var{g}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})\n\
215 Solve the set of differential-algebraic equations\n\
216 @tex\n\
217 $$ 0 = f (x, \\dot{x}, t) $$\n\
218 with\n\
219 $$ x(t_0) = x_0, \\dot{x}(t_0) = \\dot{x}_0 $$\n\
220 @end tex\n\
221 @ifnottex\n\
222 \n\
223 @example\n\
224 0 = f (x, xdot, t)\n\
225 @end example\n\
226 \n\
227 @noindent\n\
228 with\n\
229 \n\
230 @example\n\
231 x(t_0) = x_0, xdot(t_0) = xdot_0\n\
232 @end example\n\
233 \n\
234 @end ifnottex\n\
235 with functional stopping criteria (root solving).\n\
236 \n\
237 The solution is returned in the matrices @var{x} and @var{xdot},\n\
238 with each row in the result matrices corresponding to one of the\n\
239 elements in the vector @var{t_out}. The first element of @var{t}\n\
240 should be @math{t_0} and correspond to the initial state of the\n\
241 system @var{x_0} and its derivative @var{xdot_0}, so that the first\n\
242 row of the output @var{x} is @var{x_0} and the first row\n\
243 of the output @var{xdot} is @var{xdot_0}.\n\
244 \n\
245 The vector @var{t} provides an upper limit on the length of the\n\
246 integration. If the stopping condition is met, the vector\n\
247 @var{t_out} will be shorter than @var{t}, and the final element of\n\
248 @var{t_out} will be the point at which the stopping condition was met,\n\
249 and may not correspond to any element of the vector @var{t}.\n\
250 \n\
251 The first argument, @var{fcn}, is a string, inline, or function handle\n\
252 that names the function @math{f} to call to compute the vector of\n\
253 residuals for the set of equations. It must have the form\n\
254 \n\
255 @example\n\
256 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
257 @end example\n\
258 \n\
259 @noindent\n\
260 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
261 scalar.\n\
262 \n\
263 If @var{fcn} is a two-element string array or a two-element cell array\n\
264 of strings, inline functions, or function handles, the first element names\n\
265 the function @math{f} described above, and the second element names a\n\
266 function to compute the modified Jacobian\n\
267 \n\
268 @tex\n\
269 $$\n\
270 J = {\\partial f \\over \\partial x}\n\
271  + c {\\partial f \\over \\partial \\dot{x}}\n\
272 $$\n\
273 @end tex\n\
274 @ifnottex\n\
275 \n\
276 @example\n\
277 @group\n\
278  df df\n\
279 jac = -- + c ------\n\
280  dx d xdot\n\
281 @end group\n\
282 @end example\n\
283 \n\
284 @end ifnottex\n\
285 \n\
286 The modified Jacobian function must have the form\n\
287 \n\
288 @example\n\
289 @group\n\
290 \n\
291 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\
292 \n\
293 @end group\n\
294 @end example\n\
295 \n\
296 The optional second argument names a function that defines the\n\
297 constraint functions whose roots are desired during the integration.\n\
298 This function must have the form\n\
299 \n\
300 @example\n\
301 @var{g_out} = g (@var{x}, @var{t})\n\
302 @end example\n\
303 \n\
304 @noindent\n\
305 and return a vector of the constraint function values.\n\
306 If the value of any of the constraint functions changes sign, @sc{dasrt}\n\
307 will attempt to stop the integration at the point of the sign change.\n\
308 \n\
309 If the name of the constraint function is omitted, @code{dasrt} solves\n\
310 the same problem as @code{daspk} or @code{dassl}.\n\
311 \n\
312 Note that because of numerical errors in the constraint functions\n\
313 due to round-off and integration error, @sc{dasrt} may return false\n\
314 roots, or return the same root at two or more nearly equal values of\n\
315 @var{T}. If such false roots are suspected, the user should consider\n\
316 smaller error tolerances or higher precision in the evaluation of the\n\
317 constraint functions.\n\
318 \n\
319 If a root of some constraint function defines the end of the problem,\n\
320 the input to @sc{dasrt} should nevertheless allow integration to a\n\
321 point slightly past that root, so that @sc{dasrt} can locate the root\n\
322 by interpolation.\n\
323 \n\
324 The third and fourth arguments to @code{dasrt} specify the initial\n\
325 condition of the states and their derivatives, and the fourth argument\n\
326 specifies a vector of output times at which the solution is desired,\n\
327 including the time corresponding to the initial condition.\n\
328 \n\
329 The set of initial states and derivatives are not strictly required to\n\
330 be consistent. In practice, however, @sc{dassl} is not very good at\n\
331 determining a consistent set for you, so it is best if you ensure that\n\
332 the initial values result in the function evaluating to zero.\n\
333 \n\
334 The sixth argument is optional, and may be used to specify a set of\n\
335 times that the DAE solver should not integrate past. It is useful for\n\
336 avoiding difficulties with singularities and points where there is a\n\
337 discontinuity in the derivative.\n\
338 \n\
339 After a successful computation, the value of @var{istate} will be\n\
340 greater than zero (consistent with the Fortran version of @sc{dassl}).\n\
341 \n\
342 If the computation is not successful, the value of @var{istate} will be\n\
343 less than zero and @var{msg} will contain additional information.\n\
344 \n\
345 You can use the function @code{dasrt_options} to set optional\n\
346 parameters for @code{dasrt}.\n\
347 @seealso{dasrt_options, daspk, dasrt, lsode}\n\
348 @end deftypefn")
349 {
350  octave_value_list retval;
351 
352  warned_fcn_imaginary = false;
353  warned_jac_imaginary = false;
354  warned_cf_imaginary = false;
355 
356  unwind_protect frame;
357 
358  frame.protect_var (call_depth);
359  call_depth++;
360 
361  if (call_depth > 1)
362  DASRT_ABORT1 ("invalid recursive call");
363 
364  int argp = 0;
365 
366  int nargin = args.length ();
367 
368  if (nargin < 4 || nargin > 6)
369  {
370  print_usage ();
371  return retval;
372  }
373 
374  std::string fcn_name, fname, jac_name, jname;
375  dasrt_f = 0;
376  dasrt_j = 0;
377  dasrt_cf = 0;
378 
379  // Check all the arguments. Are they the right animals?
380 
381  // Here's where I take care of f and j in one shot:
382 
383  octave_value f_arg = args(0);
384 
385  if (f_arg.is_cell ())
386  {
387  Cell c = f_arg.cell_value ();
388  if (c.length () == 1)
389  f_arg = c(0);
390  else if (c.length () == 2)
391  {
392  if (c(0).is_function_handle () || c(0).is_inline_function ())
393  dasrt_f = c(0).function_value ();
394  else
395  {
396  fcn_name = unique_symbol_name ("__dasrt_fcn__");
397  fname = "function y = ";
398  fname.append (fcn_name);
399  fname.append (" (x, xdot, t) y = ");
400  dasrt_f = extract_function (c(0), "dasrt", fcn_name, fname,
401  "; endfunction");
402  }
403 
404  if (dasrt_f)
405  {
406  if (c(1).is_function_handle () || c(1).is_inline_function ())
407  dasrt_j = c(1).function_value ();
408  else
409  {
410  jac_name = unique_symbol_name ("__dasrt_jac__");
411  jname = "function jac = ";
412  jname.append (jac_name);
413  jname.append (" (x, xdot, t, cj) jac = ");
414  dasrt_j = extract_function (c(1), "dasrt", jac_name, jname,
415  "; endfunction");
416 
417  if (!dasrt_j)
418  {
419  if (fcn_name.length ())
420  clear_function (fcn_name);
421  dasrt_f = 0;
422  }
423  }
424  }
425  }
426  else
427  DASRT_ABORT1 ("incorrect number of elements in cell array");
428  }
429 
430  if (!dasrt_f && ! f_arg.is_cell ())
431  {
432  if (f_arg.is_function_handle () || f_arg.is_inline_function ())
433  dasrt_f = f_arg.function_value ();
434  else
435  {
436  switch (f_arg.rows ())
437  {
438  case 1:
439  fcn_name = unique_symbol_name ("__dasrt_fcn__");
440  fname = "function y = ";
441  fname.append (fcn_name);
442  fname.append (" (x, xdot, t) y = ");
443  dasrt_f = extract_function (f_arg, "dasrt", fcn_name, fname,
444  "; endfunction");
445  break;
446 
447  case 2:
448  {
449  string_vector tmp = args(0).all_strings ();
450 
451  if (! error_state)
452  {
453  fcn_name = unique_symbol_name ("__dasrt_fcn__");
454  fname = "function y = ";
455  fname.append (fcn_name);
456  fname.append (" (x, xdot, t) y = ");
457  dasrt_f = extract_function (tmp(0), "dasrt", fcn_name,
458  fname, "; endfunction");
459 
460  if (dasrt_f)
461  {
462  jac_name = unique_symbol_name ("__dasrt_jac__");
463  jname = "function jac = ";
464  jname.append (jac_name);
465  jname.append (" (x, xdot, t, cj) jac = ");
466  dasrt_j = extract_function (tmp(1), "dasrt", jac_name,
467  jname, "; endfunction");
468 
469  if (! dasrt_j)
470  dasrt_f = 0;
471  }
472  }
473  }
474  break;
475 
476  default:
478  ("first arg should be a string or 2-element string array");
479  }
480  }
481  }
482 
483  if (error_state || (! dasrt_f))
484  DASRT_ABORT;
485 
486  DAERTFunc func (dasrt_user_f);
487 
488  argp++;
489 
490  if (args(1).is_function_handle () || args(1).is_inline_function ())
491  {
492  dasrt_cf = args(1).function_value ();
493 
494  if (! dasrt_cf)
495  DASRT_ABORT1 ("expecting function name as argument 2");
496 
497  argp++;
498 
500  }
501  else if (args(1).is_string ())
502  {
503  dasrt_cf = is_valid_function (args(1), "dasrt", true);
504  if (! dasrt_cf)
505  DASRT_ABORT1 ("expecting function name as argument 2");
506 
507  argp++;
508 
510  }
511 
512  ColumnVector state (args(argp++).vector_value ());
513 
514  if (error_state)
515  DASRT_ABORT2 ("expecting state vector as argument %d", argp);
516 
517  ColumnVector stateprime (args(argp++).vector_value ());
518 
519  if (error_state)
520  DASRT_ABORT2 ("expecting time derivative of state vector as argument %d",
521  argp);
522 
523  ColumnVector out_times (args(argp++).vector_value ());
524 
525  if (error_state)
526  DASRT_ABORT2 ("expecting output time vector as %s argument %d", argp);
527 
528  double tzero = out_times (0);
529 
530  ColumnVector crit_times;
531 
532  bool crit_times_set = false;
533 
534  if (argp < nargin)
535  {
536  crit_times = ColumnVector (args(argp++).vector_value ());
537 
538  if (error_state)
539  DASRT_ABORT2 ("expecting critical time vector as argument %d", argp);
540 
541  crit_times_set = true;
542  }
543 
544  if (dasrt_j)
546 
547  DASRT_result output;
548 
549  DASRT dae = DASRT (state, stateprime, tzero, func);
550 
551  dae.set_options (dasrt_opts);
552 
553  if (crit_times_set)
554  output = dae.integrate (out_times, crit_times);
555  else
556  output = dae.integrate (out_times);
557 
558  if (fcn_name.length ())
559  clear_function (fcn_name);
560  if (jac_name.length ())
561  clear_function (jac_name);
562 
563  if (! error_state)
564  {
565  std::string msg = dae.error_message ();
566 
567  retval(4) = msg;
568  retval(3) = static_cast<double> (dae.integration_state ());
569 
570  if (dae.integration_ok ())
571  {
572  retval(2) = output.times ();
573  retval(1) = output.deriv ();
574  retval(0) = output.state ();
575  }
576  else
577  {
578  retval(2) = Matrix ();
579  retval(1) = Matrix ();
580  retval(0) = Matrix ();
581 
582  if (nargout < 4)
583  error ("dasrt: %s", msg.c_str ());
584  }
585  }
586 
587  return retval;
588 }
Definition: Cell.h:35
Definition: DASRT.h:69
Matrix state(void) const
Definition: DASRT.h:58
virtual octave_value_list do_multi_index_op(int nargout, const octave_value_list &idx)
Definition: ov-base.cc:203
octave_idx_type rows(void) const
Definition: ov.h:473
bool integration_ok(void) const
Definition: base-de.h:96
static uint32_t state[624]
Definition: randmtzig.c:188
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
octave_idx_type length(void) const
Definition: oct-obj.h:89
DAERTFunc & set_constraint_function(DAERTConstrFunc cf)
Definition: DAERTFunc.h:68
void protect_var(T &var)
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
static DASRT_options dasrt_opts
Definition: DASRT-opts.cc:20
static octave_function * dasrt_j
Definition: dasrt.cc:49
bool is_cell(void) const
Definition: ov.h:529
std::string error_message(void) const
Definition: DASRT.cc:566
Matrix deriv(void) const
Definition: DASRT.h:59
bool is_function_handle(void) const
Definition: ov.h:686
std::string unique_symbol_name(const std::string &basename)
Definition: variables.cc:495
Cell cell_value(void) const
Definition: ov.cc:1566
void clear_function(const std::string &nm)
Definition: variables.cc:77
#define DASRT_ABORT2(fmt, arg)
Definition: dasrt.cc:201
static octave_function * dasrt_cf
Definition: dasrt.cc:50
int error_state
Definition: error.cc:101
bool is_inline_function(void) const
Definition: ov.h:692
octave_idx_type capacity(void) const
Number of elements in the array.
Definition: Array.h:256
Definition: dMatrix.h:35
ColumnVector times(void) const
Definition: DASRT.h:60
DASRT_result integrate(const ColumnVector &tout)
Definition: DASRT.cc:388
static bool warned_cf_imaginary
Definition: dasrt.cc:55
octave_function * function_value(bool silent=false) const
Definition: ov.cc:1597
void warning(const char *fmt,...)
Definition: error.cc:681
static ColumnVector dasrt_user_f(const ColumnVector &x, const ColumnVector &xdot, double t, octave_idx_type &)
Definition: dasrt.cc:61
static ColumnVector dasrt_user_cf(const ColumnVector &x, double t)
Definition: dasrt.cc:105
static Matrix dasrt_user_j(const ColumnVector &x, const ColumnVector &xdot, double t, double cj)
Definition: dasrt.cc:145
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
static bool warned_fcn_imaginary
Definition: dasrt.cc:53
#define DASRT_ABORT
Definition: dasrt.cc:190
octave_function * extract_function(const octave_value &arg, const std::string &warn_for, const std::string &fname, const std::string &header, const std::string &trailer)
Definition: variables.cc:140
octave_idx_type integration_state(void) const
Definition: base-de.h:98
virtual octave_function * function_value(bool silent=false)
Definition: ov-base.cc:992
octave_function * is_valid_function(const std::string &fcn_name, const std::string &warn_for, bool warn)
Definition: variables.cc:99
void set_options(const DASRT_options &opt)
Definition: DASRT-opts.h:71
#define DASRT_ABORT1(msg)
Definition: dasrt.cc:193
void gripe_user_supplied_eval(const char *name)
Definition: gripes.cc:87
static octave_function * dasrt_f
Definition: dasrt.cc:48
static int call_depth
Definition: dasrt.cc:58
F77_RET_T const double * x
static bool warned_jac_imaginary
Definition: dasrt.cc:54
DAEFunc & set_jacobian_function(DAEJacFunc j)
Definition: DAEFunc.h:82