GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
daspk.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2018 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
9 the Free Software Foundation, either version 3 of the License, or
10 (at your option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License 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 <https://www.gnu.org/licenses/>.
20 
21 */
22 
23 #if defined (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 "errwarn.h"
37 #include "ovl.h"
38 #include "ov-fcn.h"
39 #include "ov-cell.h"
40 #include "pager.h"
41 #include "parse.h"
42 #include "unwind-prot.h"
43 #include "utils.h"
44 #include "variables.h"
45 
46 #include "DASPK-opts.cc"
47 
48 // Global pointer for user defined function required by daspk.
50 
51 // Global pointer for optional user defined jacobian function.
53 
54 // Have we warned about imaginary values returned from user function?
55 static bool warned_fcn_imaginary = false;
56 static bool warned_jac_imaginary = false;
57 
58 // Is this a recursive call?
59 static int call_depth = 0;
60 
63  double t, octave_idx_type& ires)
64 {
66 
67  assert (x.numel () == xdot.numel ());
68 
69  octave_value_list args;
70 
71  args(2) = t;
72  args(1) = xdot;
73  args(0) = x;
74 
75  if (daspk_fcn)
76  {
78 
79  try
80  {
81  tmp = octave::feval (daspk_fcn, args, 1);
82  }
83  catch (octave::execution_exception& e)
84  {
85  err_user_supplied_eval (e, "daspk");
86  }
87 
88  int tlen = tmp.length ();
89  if (tlen == 0 || ! tmp(0).is_defined ())
90  err_user_supplied_eval ("daspk");
91 
92  if (! warned_fcn_imaginary && tmp(0).iscomplex ())
93  {
94  warning ("daspk: ignoring imaginary part returned from user-supplied function");
95  warned_fcn_imaginary = true;
96  }
97 
98  retval = tmp(0).vector_value ();
99 
100  if (tlen > 1)
101  ires = tmp(1).idx_type_value ();
102 
103  if (retval.isempty ())
104  err_user_supplied_eval ("daspk");
105  }
106 
107  return retval;
108 }
109 
110 Matrix
112  double t, double cj)
113 {
114  Matrix retval;
115 
116  assert (x.numel () == xdot.numel ());
117 
118  octave_value_list args;
119 
120  args(3) = cj;
121  args(2) = t;
122  args(1) = xdot;
123  args(0) = x;
124 
125  if (daspk_jac)
126  {
128 
129  try
130  {
131  tmp = octave::feval (daspk_jac, args, 1);
132  }
133  catch (octave::execution_exception& e)
134  {
135  err_user_supplied_eval (e, "daspk");
136  }
137 
138  int tlen = tmp.length ();
139  if (tlen == 0 || ! tmp(0).is_defined ())
140  err_user_supplied_eval ("daspk");
141 
142  if (! warned_jac_imaginary && tmp(0).iscomplex ())
143  {
144  warning ("daspk: ignoring imaginary part returned from user-supplied jacobian function");
145  warned_jac_imaginary = true;
146  }
147 
148  retval = tmp(0).matrix_value ();
149 
150  if (retval.isempty ())
151  err_user_supplied_eval ("daspk");
152  }
153 
154  return retval;
155 }
156 
157 DEFMETHOD (daspk, interp, args, nargout,
158  doc: /* -*- texinfo -*-
159 @deftypefn {} {[@var{x}, @var{xdot}, @var{istate}, @var{msg}] =} daspk (@var{fcn}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})
160 Solve a set of differential-algebraic equations.
161 
162 @code{daspk} solves the set of equations
163 @tex
164 $$ 0 = f (x, \dot{x}, t) $$
165 with
166 $$ x(t_0) = x_0, \dot{x}(t_0) = \dot{x}_0 $$
167 @end tex
168 @ifnottex
169 
170 @example
171 0 = f (x, xdot, t)
172 @end example
173 
174 @noindent
175 with
176 
177 @example
178 x(t_0) = x_0, xdot(t_0) = xdot_0
179 @end example
180 
181 @end ifnottex
182 The solution is returned in the matrices @var{x} and @var{xdot},
183 with each row in the result matrices corresponding to one of the
184 elements in the vector @var{t}. The first element of @var{t}
185 should be @math{t_0} and correspond to the initial state of the
186 system @var{x_0} and its derivative @var{xdot_0}, so that the first
187 row of the output @var{x} is @var{x_0} and the first row
188 of the output @var{xdot} is @var{xdot_0}.
189 
190 The first argument, @var{fcn}, is a string, inline, or function handle
191 that names the function @math{f} to call to compute the vector of
192 residuals for the set of equations. It must have the form
193 
194 @example
195 @var{res} = f (@var{x}, @var{xdot}, @var{t})
196 @end example
197 
198 @noindent
199 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a
200 scalar.
201 
202 If @var{fcn} is a two-element string array or a two-element cell array
203 of strings, inline functions, or function handles, the first element names
204 the function @math{f} described above, and the second element names a
205 function to compute the modified Jacobian
206 @tex
207 $$
208 J = {\partial f \over \partial x}
209  + c {\partial f \over \partial \dot{x}}
210 $$
211 @end tex
212 @ifnottex
213 
214 @example
215 @group
216  df df
217 jac = -- + c ------
218  dx d xdot
219 @end group
220 @end example
221 
222 @end ifnottex
223 
224 The modified Jacobian function must have the form
225 
226 @example
227 @group
228 
229 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})
230 
231 @end group
232 @end example
233 
234 The second and third arguments to @code{daspk} specify the initial
235 condition of the states and their derivatives, and the fourth argument
236 specifies a vector of output times at which the solution is desired,
237 including the time corresponding to the initial condition.
238 
239 The set of initial states and derivatives are not strictly required to
240 be consistent. If they are not consistent, you must use the
241 @code{daspk_options} function to provide additional information so
242 that @code{daspk} can compute a consistent starting point.
243 
244 The fifth argument is optional, and may be used to specify a set of
245 times that the DAE solver should not integrate past. It is useful for
246 avoiding difficulties with singularities and points where there is a
247 discontinuity in the derivative.
248 
249 After a successful computation, the value of @var{istate} will be
250 greater than zero (consistent with the Fortran version of @sc{daspk}).
251 
252 If the computation is not successful, the value of @var{istate} will be
253 less than zero and @var{msg} will contain additional information.
254 
255 You can use the function @code{daspk_options} to set optional
256 parameters for @code{daspk}.
257 @seealso{dassl}
258 @end deftypefn */)
259 {
260  int nargin = args.length ();
261 
263  print_usage ();
264 
265  warned_fcn_imaginary = false;
266  warned_jac_imaginary = false;
267 
269 
271 
273  call_depth++;
274 
275  octave::symbol_table& symtab = interp.get_symbol_table ();
276 
277  if (call_depth > 1)
278  error ("daspk: invalid recursive call");
279 
280  std::string fcn_name, fname, jac_name, jname;
281  daspk_fcn = nullptr;
282  daspk_jac = nullptr;
283 
284  octave_value f_arg = args(0);
285 
286  if (f_arg.iscell ())
287  {
288  Cell c = f_arg.cell_value ();
289  if (c.numel () == 1)
290  f_arg = c(0);
291  else if (c.numel () == 2)
292  {
293  if (c(0).is_function_handle () || c(0).is_inline_function ())
294  daspk_fcn = c(0).function_value ();
295  else
296  {
297  fcn_name = unique_symbol_name ("__daspk_fcn__");
298  fname = "function y = ";
299  fname.append (fcn_name);
300  fname.append (" (x, xdot, t) y = ");
301  daspk_fcn = extract_function (c(0), "daspk", fcn_name,
302  fname, "; endfunction");
303  }
304 
305  if (daspk_fcn)
306  {
307  if (c(1).is_function_handle () || c(1).is_inline_function ())
308  daspk_jac = c(1).function_value ();
309  else
310  {
311  jac_name = unique_symbol_name ("__daspk_jac__");
312  jname = "function jac = ";
313  jname.append (jac_name);
314  jname.append (" (x, xdot, t, cj) jac = ");
315  daspk_jac = extract_function (c(1), "daspk", jac_name,
316  jname, "; endfunction");
317 
318  if (! daspk_jac)
319  {
320  if (fcn_name.length ())
321  symtab.clear_function (fcn_name);
322  daspk_fcn = nullptr;
323  }
324  }
325  }
326  }
327  else
328  error ("daspk: incorrect number of elements in cell array");
329  }
330 
331  if (! daspk_fcn && ! f_arg.iscell ())
332  {
333  if (f_arg.is_function_handle () || f_arg.is_inline_function ())
334  daspk_fcn = f_arg.function_value ();
335  else
336  {
337  switch (f_arg.rows ())
338  {
339  case 1:
340  do
341  {
342  fcn_name = unique_symbol_name ("__daspk_fcn__");
343  fname = "function y = ";
344  fname.append (fcn_name);
345  fname.append (" (x, xdot, t) y = ");
346  daspk_fcn = extract_function (f_arg, "daspk", fcn_name,
347  fname, "; endfunction");
348  }
349  while (0);
350  break;
351 
352  case 2:
353  {
355 
356  fcn_name = unique_symbol_name ("__daspk_fcn__");
357  fname = "function y = ";
358  fname.append (fcn_name);
359  fname.append (" (x, xdot, t) y = ");
360  daspk_fcn = extract_function (tmp(0), "daspk", fcn_name,
361  fname, "; endfunction");
362 
363  if (daspk_fcn)
364  {
365  jac_name = unique_symbol_name ("__daspk_jac__");
366  jname = "function jac = ";
367  jname.append (jac_name);
368  jname.append (" (x, xdot, t, cj) jac = ");
369  daspk_jac = extract_function (tmp(1), "daspk", jac_name,
370  jname, "; endfunction");
371 
372  if (! daspk_jac)
373  {
374  if (fcn_name.length ())
375  symtab.clear_function (fcn_name);
376  daspk_fcn = nullptr;
377  }
378  }
379  }
380  }
381  }
382  }
383 
384  if (! daspk_fcn)
385  return retval;
386 
387  ColumnVector state = args(1).xvector_value ("daspk: initial state X_0 must be a vector");
388 
389  ColumnVector deriv = args(2).xvector_value ("daspk: initial derivatives XDOT_0 must be a vector");
390 
391  ColumnVector out_times = args(3).xvector_value ("daspk: output time variable T must be a vector");
392 
393  ColumnVector crit_times;
394  int crit_times_set = 0;
395  if (nargin > 4)
396  {
397  crit_times = args(4).xvector_value ("daspk: list of critical times T_CRIT must be a vector");
398 
399  crit_times_set = 1;
400  }
401 
402  if (state.numel () != deriv.numel ())
403  error ("daspk: X_0 and XDOT_0 must have the same size");
404 
405  double tzero = out_times (0);
406 
408  if (daspk_jac)
410 
411  DASPK dae (state, deriv, tzero, func);
412  dae.set_options (daspk_opts);
413 
414  Matrix output;
415  Matrix deriv_output;
416 
417  if (crit_times_set)
418  output = dae.integrate (out_times, deriv_output, crit_times);
419  else
420  output = dae.integrate (out_times, deriv_output);
421 
422  if (fcn_name.length ())
423  symtab.clear_function (fcn_name);
424  if (jac_name.length ())
425  symtab.clear_function (jac_name);
426 
427  std::string msg = dae.error_message ();
428 
429  if (dae.integration_ok ())
430  {
431  retval(0) = output;
432  retval(1) = deriv_output;
433  }
434  else
435  {
436  if (nargout < 3)
437  error ("daspk: %s", msg.c_str ());
438 
439  retval(0) = Matrix ();
440  retval(1) = Matrix ();
441  }
442 
443  retval(2) = static_cast<double> (dae.integration_state ());
444  retval(3) = msg;
445 
446  return retval;
447 }
OCTINTERP_API octave_value_list feval(const std::string &name, const octave_value_list &args=octave_value_list(), int nargout=0)
Definition: Cell.h:37
#define DEFMETHOD(name, interp_name, args_name, nargout_name, doc)
Macro to define a builtin method.
Definition: defun.h:135
Matrix integrate(const ColumnVector &tout, Matrix &xdot_out)
Definition: DASPK.cc:551
fname
Definition: load-save.cc:767
bool isempty(void) const
Definition: ov.h:529
OCTINTERP_API void print_usage(void)
Definition: defun.cc:54
std::string error_message(void) const
Definition: DASPK.cc:699
static octave_function * daspk_fcn
Definition: daspk.cc:49
void error(const char *fmt,...)
Definition: error.cc:578
octave_idx_type integration_state(void) const
Definition: base-de.h:99
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function t
Definition: ov-usr-fcn.cc:997
nd example oindent opens the file binary numeric values will be read assuming they are stored in IEEE format with the least significant bit and then converted to the native representation Opening a file that is already open simply opens it again and returns a separate file id It is not an error to open a file several though writing to the same file through several different file ids may produce unexpected results The possible values of text mode reading and writing automatically converts linefeeds to the appropriate line end character for the you may append a you must also open the file in binary mode The parameter conversions are currently only supported for and permissions will be set to and then everything is written in a single operation This is very efficient and improves performance c
Definition: file-io.cc:587
i e
Definition: data.cc:2591
static bool warned_fcn_imaginary
Definition: daspk.cc:55
std::string unique_symbol_name(const std::string &basename)
Definition: variables.cc:498
bool integration_ok(void) const
Definition: base-de.h:97
Matrix daspk_user_jacobian(const ColumnVector &x, const ColumnVector &xdot, double t, double cj)
Definition: daspk.cc:111
Definition: DASPK.h:35
OCTAVE_EXPORT octave_value_list return the number of command line arguments passed to Octave If called with the optional argument the function xample nargout(@histc)
Definition: ov-usr-fcn.cc:997
bool is_function_handle(void) const
Definition: ov.h:749
bool iscell(void) const
Definition: ov.h:536
octave_idx_type rows(void) const
Definition: ov.h:472
double tmp
Definition: data.cc:6252
octave_value retval
Definition: data.cc:6246
octave_function * function_value(bool silent=false) const
Definition: dMatrix.h:36
static uint32_t state[624]
Definition: randmtzig.cc:183
void warning(const char *fmt,...)
Definition: error.cc:801
octave::unwind_protect frame
Definition: graphics.cc:12190
void err_user_supplied_eval(const char *name)
Definition: errwarn.cc:148
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:125
static bool warned_jac_imaginary
Definition: daspk.cc:56
args.length() nargin
Definition: file-io.cc:589
ColumnVector daspk_user_function(const ColumnVector &x, const ColumnVector &xdot, double t, octave_idx_type &ires)
Definition: daspk.cc:62
void clear_function(const std::string &name)
Definition: symtab.h:392
string_vector string_vector_value(bool pad=false) const
Definition: ov.h:958
static octave_function * daspk_jac
Definition: daspk.cc:52
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:366
static int call_depth
Definition: daspk.cc:59
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:888
Cell cell_value(void) const
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x
bool is_inline_function(void) const
Definition: ov.h:755
DAEFunc & set_jacobian_function(DAEJacFunc j)
Definition: DAEFunc.h:84