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
daspk.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 "DASPK.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 "DASPK-opts.cc"
46 
47 // Global pointer for user defined function required by daspk.
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 (daspk_fcn)
75  {
76  octave_value_list tmp = daspk_fcn->do_multi_index_op (1, args);
77 
78  if (error_state)
79  {
80  gripe_user_supplied_eval ("daspk");
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 ("daspk: 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 ("daspk");
100  }
101  else
102  gripe_user_supplied_eval ("daspk");
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 (daspk_jac)
124  {
125  octave_value_list tmp = daspk_jac->do_multi_index_op (1, args);
126 
127  if (error_state)
128  {
129  gripe_user_supplied_eval ("daspk");
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 ("daspk: 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 ("daspk");
146  }
147  else
148  gripe_user_supplied_eval ("daspk");
149  }
150 
151  return retval;
152 }
153 
154 #define DASPK_ABORT() \
155  return retval
156 
157 #define DASPK_ABORT1(msg) \
158  do \
159  { \
160  ::error ("daspk: " msg); \
161  DASPK_ABORT (); \
162  } \
163  while (0)
164 
165 #define DASPK_ABORT2(fmt, arg) \
166  do \
167  { \
168  ::error ("daspk: " fmt, arg); \
169  DASPK_ABORT (); \
170  } \
171  while (0)
172 
173 DEFUN (daspk, args, nargout,
174  "-*- texinfo -*-\n\
175 @deftypefn {Built-in Function} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} daspk (@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 @tex\n\
221 $$\n\
222 J = {\\partial f \\over \\partial x}\n\
223  + c {\\partial f \\over \\partial \\dot{x}}\n\
224 $$\n\
225 @end tex\n\
226 @ifnottex\n\
227 \n\
228 @example\n\
229 @group\n\
230  df df\n\
231 jac = -- + c ------\n\
232  dx d xdot\n\
233 @end group\n\
234 @end example\n\
235 \n\
236 @end ifnottex\n\
237 \n\
238 The modified Jacobian function must have the form\n\
239 \n\
240 @example\n\
241 @group\n\
242 \n\
243 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\
244 \n\
245 @end group\n\
246 @end example\n\
247 \n\
248 The second and third arguments to @code{daspk} specify the initial\n\
249 condition of the states and their derivatives, and the fourth argument\n\
250 specifies a vector of output times at which the solution is desired,\n\
251 including the time corresponding to the initial condition.\n\
252 \n\
253 The set of initial states and derivatives are not strictly required to\n\
254 be consistent. If they are not consistent, you must use the\n\
255 @code{daspk_options} function to provide additional information so\n\
256 that @code{daspk} can compute a consistent starting point.\n\
257 \n\
258 The fifth argument is optional, and may be used to specify a set of\n\
259 times that the DAE solver should not integrate past. It is useful for\n\
260 avoiding difficulties with singularities and points where there is a\n\
261 discontinuity in the derivative.\n\
262 \n\
263 After a successful computation, the value of @var{istate} will be\n\
264 greater than zero (consistent with the Fortran version of @sc{daspk}).\n\
265 \n\
266 If the computation is not successful, the value of @var{istate} will be\n\
267 less than zero and @var{msg} will contain additional information.\n\
268 \n\
269 You can use the function @code{daspk_options} to set optional\n\
270 parameters for @code{daspk}.\n\
271 @seealso{dassl}\n\
272 @end deftypefn")
273 {
274  octave_value_list retval;
275 
276  warned_fcn_imaginary = false;
277  warned_jac_imaginary = false;
278 
279  unwind_protect frame;
280 
281  frame.protect_var (call_depth);
282  call_depth++;
283 
284  if (call_depth > 1)
285  DASPK_ABORT1 ("invalid recursive call");
286 
287  int nargin = args.length ();
288 
289  if (nargin > 3 && nargin < 6)
290  {
291  std::string fcn_name, fname, jac_name, jname;
292  daspk_fcn = 0;
293  daspk_jac = 0;
294 
295  octave_value f_arg = args(0);
296 
297  if (f_arg.is_cell ())
298  {
299  Cell c = f_arg.cell_value ();
300  if (c.length () == 1)
301  f_arg = c(0);
302  else if (c.length () == 2)
303  {
304  if (c(0).is_function_handle () || c(0).is_inline_function ())
305  daspk_fcn = c(0).function_value ();
306  else
307  {
308  fcn_name = unique_symbol_name ("__daspk_fcn__");
309  fname = "function y = ";
310  fname.append (fcn_name);
311  fname.append (" (x, xdot, t) y = ");
312  daspk_fcn = extract_function
313  (c(0), "daspk", fcn_name, fname, "; endfunction");
314  }
315 
316  if (daspk_fcn)
317  {
318  if (c(1).is_function_handle () || c(1).is_inline_function ())
319  daspk_jac = c(1).function_value ();
320  else
321  {
322  jac_name = unique_symbol_name ("__daspk_jac__");
323  jname = "function jac = ";
324  jname.append (jac_name);
325  jname.append (" (x, xdot, t, cj) jac = ");
326  daspk_jac = extract_function (c(1), "daspk", jac_name,
327  jname, "; endfunction");
328 
329  if (!daspk_jac)
330  {
331  if (fcn_name.length ())
332  clear_function (fcn_name);
333  daspk_fcn = 0;
334  }
335  }
336  }
337  }
338  else
339  DASPK_ABORT1 ("incorrect number of elements in cell array");
340  }
341 
342  if (!daspk_fcn && ! f_arg.is_cell ())
343  {
344  if (f_arg.is_function_handle () || f_arg.is_inline_function ())
345  daspk_fcn = f_arg.function_value ();
346  else
347  {
348  switch (f_arg.rows ())
349  {
350  case 1:
351  do
352  {
353  fcn_name = unique_symbol_name ("__daspk_fcn__");
354  fname = "function y = ";
355  fname.append (fcn_name);
356  fname.append (" (x, xdot, t) y = ");
357  daspk_fcn = extract_function (f_arg, "daspk", fcn_name,
358  fname, "; endfunction");
359  }
360  while (0);
361  break;
362 
363  case 2:
364  {
365  string_vector tmp = f_arg.all_strings ();
366 
367  if (! error_state)
368  {
369  fcn_name = unique_symbol_name ("__daspk_fcn__");
370  fname = "function y = ";
371  fname.append (fcn_name);
372  fname.append (" (x, xdot, t) y = ");
373  daspk_fcn = extract_function (tmp(0), "daspk", fcn_name,
374  fname, "; endfunction");
375 
376  if (daspk_fcn)
377  {
378  jac_name = unique_symbol_name ("__daspk_jac__");
379  jname = "function jac = ";
380  jname.append (jac_name);
381  jname.append (" (x, xdot, t, cj) jac = ");
382  daspk_jac = extract_function (tmp(1), "daspk",
383  jac_name, jname,
384  "; endfunction");
385 
386  if (!daspk_jac)
387  {
388  if (fcn_name.length ())
389  clear_function (fcn_name);
390  daspk_fcn = 0;
391  }
392  }
393  }
394  }
395  }
396  }
397  }
398 
399  if (error_state || ! daspk_fcn)
400  DASPK_ABORT ();
401 
402  ColumnVector state = ColumnVector (args(1).vector_value ());
403 
404  if (error_state)
405  DASPK_ABORT1 ("expecting state vector as second argument");
406 
407  ColumnVector deriv (args(2).vector_value ());
408 
409  if (error_state)
410  DASPK_ABORT1 ("expecting derivative vector as third argument");
411 
412  ColumnVector out_times (args(3).vector_value ());
413 
414  if (error_state)
415  DASPK_ABORT1 ("expecting output time vector as fourth argument");
416 
417  ColumnVector crit_times;
418  int crit_times_set = 0;
419  if (nargin > 4)
420  {
421  crit_times = ColumnVector (args(4).vector_value ());
422 
423  if (error_state)
424  DASPK_ABORT1 ("expecting critical time vector as fifth argument");
425 
426  crit_times_set = 1;
427  }
428 
429  if (state.capacity () != deriv.capacity ())
430  DASPK_ABORT1 ("x and xdot must have the same size");
431 
432  double tzero = out_times (0);
433 
435  if (daspk_jac)
437 
438  DASPK dae (state, deriv, tzero, func);
439  dae.set_options (daspk_opts);
440 
441  Matrix output;
442  Matrix deriv_output;
443 
444  if (crit_times_set)
445  output = dae.integrate (out_times, deriv_output, crit_times);
446  else
447  output = dae.integrate (out_times, deriv_output);
448 
449  if (fcn_name.length ())
450  clear_function (fcn_name);
451  if (jac_name.length ())
452  clear_function (jac_name);
453 
454  if (! error_state)
455  {
456  std::string msg = dae.error_message ();
457 
458  retval(3) = msg;
459  retval(2) = static_cast<double> (dae.integration_state ());
460 
461  if (dae.integration_ok ())
462  {
463  retval(1) = deriv_output;
464  retval(0) = output;
465  }
466  else
467  {
468  retval(1) = Matrix ();
469  retval(0) = Matrix ();
470 
471  if (nargout < 3)
472  error ("daspk: %s", msg.c_str ());
473  }
474  }
475  }
476  else
477  print_usage ();
478 
479  return retval;
480 }