GNU Octave  3.8.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
dassl.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2013 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 <string>
28 
29 #include <iomanip>
30 #include <iostream>
31 
32 #include "DASSL.h"
33 
34 #include "defun.h"
35 #include "error.h"
36 #include "gripes.h"
37 #include "oct-obj.h"
38 #include "ov-fcn.h"
39 #include "ov-cell.h"
40 #include "pager.h"
41 #include "unwind-prot.h"
42 #include "utils.h"
43 #include "variables.h"
44 
45 #include "DASSL-opts.cc"
46 
47 // Global pointer for user defined function required by dassl.
49 
50 // Global pointer for optional user defined jacobian function.
52 
53 // Have we warned about imaginary values returned from user function?
54 static bool warned_fcn_imaginary = false;
55 static bool warned_jac_imaginary = false;
56 
57 // Is this a recursive call?
58 static int call_depth = 0;
59 
62  double t, octave_idx_type& ires)
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 (dassl_fcn)
75  {
76  octave_value_list tmp = dassl_fcn->do_multi_index_op (1, args);
77 
78  if (error_state)
79  {
80  gripe_user_supplied_eval ("dassl");
81  return retval;
82  }
83 
84  int tlen = tmp.length ();
85  if (tlen > 0 && tmp(0).is_defined ())
86  {
87  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
88  {
89  warning ("dassl: ignoring imaginary part returned from user-supplied function");
90  warned_fcn_imaginary = true;
91  }
92 
93  retval = ColumnVector (tmp(0).vector_value ());
94 
95  if (tlen > 1)
96  ires = tmp(1).int_value ();
97 
98  if (error_state || retval.length () == 0)
99  gripe_user_supplied_eval ("dassl");
100  }
101  else
102  gripe_user_supplied_eval ("dassl");
103  }
104 
105  return retval;
106 }
107 
108 Matrix
110  double t, double cj)
111 {
112  Matrix retval;
113 
114  assert (x.capacity () == xdot.capacity ());
115 
116  octave_value_list args;
117 
118  args(3) = cj;
119  args(2) = t;
120  args(1) = xdot;
121  args(0) = x;
122 
123  if (dassl_jac)
124  {
125  octave_value_list tmp = dassl_jac->do_multi_index_op (1, args);
126 
127  if (error_state)
128  {
129  gripe_user_supplied_eval ("dassl");
130  return retval;
131  }
132 
133  int tlen = tmp.length ();
134  if (tlen > 0 && tmp(0).is_defined ())
135  {
136  if (! warned_jac_imaginary && tmp(0).is_complex_type ())
137  {
138  warning ("dassl: ignoring imaginary part returned from user-supplied jacobian function");
139  warned_jac_imaginary = true;
140  }
141 
142  retval = tmp(0).matrix_value ();
143 
144  if (error_state || retval.length () == 0)
145  gripe_user_supplied_eval ("dassl");
146  }
147  else
148  gripe_user_supplied_eval ("dassl");
149  }
150 
151  return retval;
152 }
153 
154 #define DASSL_ABORT() \
155  return retval
156 
157 #define DASSL_ABORT1(msg) \
158  do \
159  { \
160  ::error ("dassl: " msg); \
161  DASSL_ABORT (); \
162  } \
163  while (0)
164 
165 #define DASSL_ABORT2(fmt, arg) \
166  do \
167  { \
168  ::error ("dassl: " fmt, arg); \
169  DASSL_ABORT (); \
170  } \
171  while (0)
172 
173 DEFUN (dassl, args, nargout,
174  "-*- texinfo -*-\n\
175 @deftypefn {Built-in Function} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} dassl (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})\n\
176 Solve the set of differential-algebraic equations\n\
177 @tex\n\
178 $$ 0 = f (x, \\dot{x}, t) $$\n\
179 with\n\
180 $$ x(t_0) = x_0, \\dot{x}(t_0) = \\dot{x}_0 $$\n\
181 @end tex\n\
182 @ifnottex\n\
183 \n\
184 @example\n\
185 0 = f (x, xdot, t)\n\
186 @end example\n\
187 \n\
188 @noindent\n\
189 with\n\
190 \n\
191 @example\n\
192 x(t_0) = x_0, xdot(t_0) = xdot_0\n\
193 @end example\n\
194 \n\
195 @end ifnottex\n\
196 The solution is returned in the matrices @var{x} and @var{xdot},\n\
197 with each row in the result matrices corresponding to one of the\n\
198 elements in the vector @var{t}. The first element of @var{t}\n\
199 should be @math{t_0} and correspond to the initial state of the\n\
200 system @var{x_0} and its derivative @var{xdot_0}, so that the first\n\
201 row of the output @var{x} is @var{x_0} and the first row\n\
202 of the output @var{xdot} is @var{xdot_0}.\n\
203 \n\
204 The first argument, @var{fcn}, is a string, inline, or function handle\n\
205 that names the function @math{f} to call to compute the vector of\n\
206 residuals for the set of equations. It must have the form\n\
207 \n\
208 @example\n\
209 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
210 @end example\n\
211 \n\
212 @noindent\n\
213 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
214 scalar.\n\
215 \n\
216 If @var{fcn} is a two-element string array or a two-element cell array\n\
217 of strings, inline functions, or function handles, the first element names\n\
218 the function @math{f} described above, and the second element names a\n\
219 function to compute the modified Jacobian\n\
220 \n\
221 @tex\n\
222 $$\n\
223 J = {\\partial f \\over \\partial x}\n\
224  + c {\\partial f \\over \\partial \\dot{x}}\n\
225 $$\n\
226 @end tex\n\
227 @ifnottex\n\
228 \n\
229 @example\n\
230 @group\n\
231  df df\n\
232 jac = -- + c ------\n\
233  dx d xdot\n\
234 @end group\n\
235 @end example\n\
236 \n\
237 @end ifnottex\n\
238 \n\
239 The modified Jacobian function must have the form\n\
240 \n\
241 @example\n\
242 @group\n\
243 \n\
244 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\
245 \n\
246 @end group\n\
247 @end example\n\
248 \n\
249 The second and third arguments to @code{dassl} specify the initial\n\
250 condition of the states and their derivatives, and the fourth argument\n\
251 specifies a vector of output times at which the solution is desired,\n\
252 including the time corresponding to the initial condition.\n\
253 \n\
254 The set of initial states and derivatives are not strictly required to\n\
255 be consistent. In practice, however, @sc{dassl} is not very good at\n\
256 determining a consistent set for you, so it is best if you ensure that\n\
257 the initial values result in the function evaluating to zero.\n\
258 \n\
259 The fifth argument is optional, and may be used to specify a set of\n\
260 times that the DAE solver should not integrate past. It is useful for\n\
261 avoiding difficulties with singularities and points where there is a\n\
262 discontinuity in the derivative.\n\
263 \n\
264 After a successful computation, the value of @var{istate} will be\n\
265 greater than zero (consistent with the Fortran version of @sc{dassl}).\n\
266 \n\
267 If the computation is not successful, the value of @var{istate} will be\n\
268 less than zero and @var{msg} will contain additional information.\n\
269 \n\
270 You can use the function @code{dassl_options} to set optional\n\
271 parameters for @code{dassl}.\n\
272 @seealso{daspk, dasrt, lsode}\n\
273 @end deftypefn")
274 {
275  octave_value_list retval;
276 
277  warned_fcn_imaginary = false;
278  warned_jac_imaginary = false;
279 
280  unwind_protect frame;
281 
282  frame.protect_var (call_depth);
283  call_depth++;
284 
285  if (call_depth > 1)
286  DASSL_ABORT1 ("invalid recursive call");
287 
288  int nargin = args.length ();
289 
290  if (nargin > 3 && nargin < 6 && nargout < 5)
291  {
292  std::string fcn_name, fname, jac_name, jname;
293  dassl_fcn = 0;
294  dassl_jac = 0;
295 
296  octave_value f_arg = args(0);
297 
298  if (f_arg.is_cell ())
299  {
300  Cell c = f_arg.cell_value ();
301  if (c.length () == 1)
302  f_arg = c(0);
303  else if (c.length () == 2)
304  {
305  if (c(0).is_function_handle () || c(0).is_inline_function ())
306  dassl_fcn = c(0).function_value ();
307  else
308  {
309  fcn_name = unique_symbol_name ("__dassl_fcn__");
310  fname = "function y = ";
311  fname.append (fcn_name);
312  fname.append (" (x, xdot, t) y = ");
313  dassl_fcn = extract_function (c(0), "dassl", fcn_name, fname,
314  "; endfunction");
315  }
316 
317  if (dassl_fcn)
318  {
319  if (c(1).is_function_handle () || c(1).is_inline_function ())
320  dassl_jac = c(1).function_value ();
321  else
322  {
323  jac_name = unique_symbol_name ("__dassl_jac__");
324  jname = "function jac = ";
325  jname.append (jac_name);
326  jname.append (" (x, xdot, t, cj) jac = ");
327  dassl_jac = extract_function (c(1), "dassl", jac_name,
328  jname, "; endfunction");
329 
330  if (!dassl_jac)
331  {
332  if (fcn_name.length ())
333  clear_function (fcn_name);
334  dassl_fcn = 0;
335  }
336  }
337  }
338  }
339  else
340  DASSL_ABORT1 ("incorrect number of elements in cell array");
341  }
342 
343  if (!dassl_fcn && ! f_arg.is_cell ())
344  {
345  if (f_arg.is_function_handle () || f_arg.is_inline_function ())
346  dassl_fcn = f_arg.function_value ();
347  else
348  {
349  switch (f_arg.rows ())
350  {
351  case 1:
352  do
353  {
354  fcn_name = unique_symbol_name ("__dassl_fcn__");
355  fname = "function y = ";
356  fname.append (fcn_name);
357  fname.append (" (x, xdot, t) y = ");
358  dassl_fcn = extract_function (f_arg, "dassl", fcn_name,
359  fname, "; endfunction");
360  }
361  while (0);
362  break;
363 
364  case 2:
365  {
366  string_vector tmp = f_arg.all_strings ();
367 
368  if (! error_state)
369  {
370  fcn_name = unique_symbol_name ("__dassl_fcn__");
371  fname = "function y = ";
372  fname.append (fcn_name);
373  fname.append (" (x, xdot, t) y = ");
374  dassl_fcn = extract_function (tmp(0), "dassl", fcn_name,
375  fname, "; endfunction");
376 
377  if (dassl_fcn)
378  {
379  jac_name = unique_symbol_name ("__dassl_jac__");
380  jname = "function jac = ";
381  jname.append (jac_name);
382  jname.append (" (x, xdot, t, cj) jac = ");
383  dassl_jac = extract_function (tmp(1), "dassl",
384  jac_name, jname,
385  "; endfunction");
386 
387  if (!dassl_jac)
388  {
389  if (fcn_name.length ())
390  clear_function (fcn_name);
391  dassl_fcn = 0;
392  }
393  }
394  }
395  }
396  }
397  }
398  }
399 
400  if (error_state || ! dassl_fcn)
401  DASSL_ABORT ();
402 
403  ColumnVector state = ColumnVector (args(1).vector_value ());
404 
405  if (error_state)
406  DASSL_ABORT1 ("expecting state vector as second argument");
407 
408  ColumnVector deriv (args(2).vector_value ());
409 
410  if (error_state)
411  DASSL_ABORT1 ("expecting derivative vector as third argument");
412 
413  ColumnVector out_times (args(3).vector_value ());
414 
415  if (error_state)
416  DASSL_ABORT1 ("expecting output time vector as fourth argument");
417 
418  ColumnVector crit_times;
419  int crit_times_set = 0;
420  if (nargin > 4)
421  {
422  crit_times = ColumnVector (args(4).vector_value ());
423 
424  if (error_state)
425  DASSL_ABORT1 ("expecting critical time vector as fifth argument");
426 
427  crit_times_set = 1;
428  }
429 
430  if (state.capacity () != deriv.capacity ())
431  DASSL_ABORT1 ("x and xdot must have the same size");
432 
433  double tzero = out_times (0);
434 
436  if (dassl_jac)
438 
439  DASSL dae (state, deriv, tzero, func);
440 
441  dae.set_options (dassl_opts);
442 
443  Matrix output;
444  Matrix deriv_output;
445 
446  if (crit_times_set)
447  output = dae.integrate (out_times, deriv_output, crit_times);
448  else
449  output = dae.integrate (out_times, deriv_output);
450 
451  if (fcn_name.length ())
452  clear_function (fcn_name);
453  if (jac_name.length ())
454  clear_function (jac_name);
455 
456  if (! error_state)
457  {
458  std::string msg = dae.error_message ();
459 
460  retval(3) = msg;
461  retval(2) = static_cast<double> (dae.integration_state ());
462 
463  if (dae.integration_ok ())
464  {
465  retval(1) = deriv_output;
466  retval(0) = output;
467  }
468  else
469  {
470  retval(1) = Matrix ();
471  retval(0) = Matrix ();
472 
473  if (nargout < 3)
474  error ("dassl: %s", msg.c_str ());
475  }
476  }
477  }
478  else
479  print_usage ();
480 
481  return retval;
482 }
483 
484 /*
485 ## dassl-1.m
486 ##
487 ## Test dassl() function
488 ##
489 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
490 ## Comalco Research and Technology
491 ## 20 May 1998
492 ##
493 ## Problem
494 ##
495 ## y1' = -y2, y1(0) = 1
496 ## y2' = y1, y2(0) = 0
497 ##
498 ## Solution
499 ##
500 ## y1(t) = cos(t)
501 ## y2(t) = sin(t)
502 ##
503 %!function res = __f (x, xdot, t)
504 %! res = [xdot(1)+x(2); xdot(2)-x(1)];
505 %!endfunction
506 
507 %!test
508 %!
509 %! x0 = [1; 0];
510 %! xdot0 = [0; 1];
511 %! t = (0:1:10)';
512 %!
513 %! tol = 100 * dassl_options ("relative tolerance");
514 %!
515 %! [x, xdot] = dassl ("__f", x0, xdot0, t);
516 %!
517 %! y = [cos(t), sin(t)];
518 %!
519 %! assert (x, y, tol);
520 
521 ## dassl-2.m
522 ##
523 ## Test dassl() function
524 ##
525 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
526 ## Comalco Research and Technology
527 ## 20 May 1998
528 ##
529 ## Based on SLATEC quick check for DASSL by Linda Petzold
530 ##
531 ## Problem
532 ##
533 ## x1' + 10*x1 = 0, x1(0) = 1
534 ## x1 + x2 = 1, x2(0) = 0
535 ##
536 ##
537 ## Solution
538 ##
539 ## x1(t) = exp(-10*t)
540 ## x2(t) = 1 - x(1)
541 ##
542 %!function res = __f (x, xdot, t)
543 %! res = [xdot(1)+10*x(1); x(1)+x(2)-1];
544 %!endfunction
545 
546 %!test
547 %!
548 %! x0 = [1; 0];
549 %! xdot0 = [-10; 10];
550 %! t = (0:0.2:1)';
551 %!
552 %! tol = 500 * dassl_options ("relative tolerance");
553 %!
554 %! [x, xdot] = dassl ("__f", x0, xdot0, t);
555 %!
556 %! y = [exp(-10*t), 1-exp(-10*t)];
557 %!
558 %! assert (x, y, tol);
559 
560 %!test
561 %! old_tol = dassl_options ("absolute tolerance");
562 %! dassl_options ("absolute tolerance", eps);
563 %! assert (dassl_options ("absolute tolerance") == eps);
564 %! ## Restore old value of tolerance
565 %! dassl_options ("absolute tolerance", old_tol);
566 
567 %!error dassl_options ("foo", 1, 2)
568 */