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
lsode.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-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 <string>
28 
29 #include <iomanip>
30 #include <iostream>
31 
32 #include "LSODE.h"
33 #include "lo-mappers.h"
34 
35 #include "defun.h"
36 #include "error.h"
37 #include "gripes.h"
38 #include "oct-obj.h"
39 #include "ov-fcn.h"
40 #include "ov-cell.h"
41 #include "pager.h"
42 #include "pr-output.h"
43 #include "unwind-prot.h"
44 #include "utils.h"
45 #include "variables.h"
46 
47 #include "LSODE-opts.cc"
48 
49 // Global pointer for user defined function required by lsode.
51 
52 // Global pointer for optional user defined jacobian function used by lsode.
54 
55 // Have we warned about imaginary values returned from user function?
56 static bool warned_fcn_imaginary = false;
57 static bool warned_jac_imaginary = false;
58 
59 // Is this a recursive call?
60 static int call_depth = 0;
61 
64 {
65  ColumnVector retval;
66 
67  octave_value_list args;
68  args(1) = t;
69  args(0) = x;
70 
71  if (lsode_fcn)
72  {
73  octave_value_list tmp = lsode_fcn->do_multi_index_op (1, args);
74 
75  if (error_state)
76  {
77  gripe_user_supplied_eval ("lsode");
78  return retval;
79  }
80 
81  if (tmp.length () > 0 && tmp(0).is_defined ())
82  {
83  if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
84  {
85  warning ("lsode: ignoring imaginary part returned from user-supplied function");
86  warned_fcn_imaginary = true;
87  }
88 
89  retval = ColumnVector (tmp(0).vector_value ());
90 
91  if (error_state || retval.length () == 0)
92  gripe_user_supplied_eval ("lsode");
93  }
94  else
95  gripe_user_supplied_eval ("lsode");
96  }
97 
98  return retval;
99 }
100 
101 Matrix
103 {
104  Matrix retval;
105 
106  octave_value_list args;
107  args(1) = t;
108  args(0) = x;
109 
110  if (lsode_jac)
111  {
112  octave_value_list tmp = lsode_jac->do_multi_index_op (1, args);
113 
114  if (error_state)
115  {
116  gripe_user_supplied_eval ("lsode");
117  return retval;
118  }
119 
120  if (tmp.length () > 0 && tmp(0).is_defined ())
121  {
122  if (! warned_jac_imaginary && tmp(0).is_complex_type ())
123  {
124  warning ("lsode: ignoring imaginary part returned from user-supplied jacobian function");
125  warned_jac_imaginary = true;
126  }
127 
128  retval = tmp(0).matrix_value ();
129 
130  if (error_state || retval.length () == 0)
131  gripe_user_supplied_eval ("lsode");
132  }
133  else
134  gripe_user_supplied_eval ("lsode");
135  }
136 
137  return retval;
138 }
139 
140 #define LSODE_ABORT() \
141  return retval
142 
143 #define LSODE_ABORT1(msg) \
144  do \
145  { \
146  ::error ("lsode: " msg); \
147  LSODE_ABORT (); \
148  } \
149  while (0)
150 
151 #define LSODE_ABORT2(fmt, arg) \
152  do \
153  { \
154  ::error ("lsode: " fmt, arg); \
155  LSODE_ABORT (); \
156  } \
157  while (0)
158 
159 DEFUN (lsode, args, nargout,
160  "-*- texinfo -*-\n\
161 @deftypefn {Built-in Function} {[@var{x}, @var{istate}, @var{msg}] =} lsode (@var{fcn}, @var{x_0}, @var{t})\n\
162 @deftypefnx {Built-in Function} {[@var{x}, @var{istate}, @var{msg}] =} lsode (@var{fcn}, @var{x_0}, @var{t}, @var{t_crit})\n\
163 Ordinary Differential Equation (ODE) solver.\n\
164 \n\
165 The set of differential equations to solve is\n\
166 @tex\n\
167 $$ {dx \\over dt} = f (x, t) $$\n\
168 with\n\
169 $$ x(t_0) = x_0 $$\n\
170 @end tex\n\
171 @ifnottex\n\
172 \n\
173 @example\n\
174 @group\n\
175 dx\n\
176 -- = f (x, t)\n\
177 dt\n\
178 @end group\n\
179 @end example\n\
180 \n\
181 @noindent\n\
182 with\n\
183 \n\
184 @example\n\
185 x(t_0) = x_0\n\
186 @end example\n\
187 \n\
188 @end ifnottex\n\
189 The solution is returned in the matrix @var{x}, with each row\n\
190 corresponding to an element of the vector @var{t}. The first element\n\
191 of @var{t} should be @math{t_0} and should correspond to the initial\n\
192 state of the system @var{x_0}, so that the first row of the output\n\
193 is @var{x_0}.\n\
194 \n\
195 The first argument, @var{fcn}, is a string, inline, or function handle\n\
196 that names the function @math{f} to call to compute the vector of right\n\
197 hand sides for the set of equations. The function must have the form\n\
198 \n\
199 @example\n\
200 @var{xdot} = f (@var{x}, @var{t})\n\
201 @end example\n\
202 \n\
203 @noindent\n\
204 in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar.\n\
205 \n\
206 If @var{fcn} is a two-element string array or a two-element cell array\n\
207 of strings, inline functions, or function handles, the first element names\n\
208 the function @math{f} described above, and the second element names a\n\
209 function to compute the Jacobian of @math{f}. The Jacobian function\n\
210 must have the form\n\
211 \n\
212 @example\n\
213 @var{jac} = j (@var{x}, @var{t})\n\
214 @end example\n\
215 \n\
216 @noindent\n\
217 in which @var{jac} is the matrix of partial derivatives\n\
218 @tex\n\
219 $$ J = {\\partial f_i \\over \\partial x_j} = \\left[\\matrix{\n\
220 {\\partial f_1 \\over \\partial x_1}\n\
221  & {\\partial f_1 \\over \\partial x_2}\n\
222  & \\cdots\n\
223  & {\\partial f_1 \\over \\partial x_N} \\cr\n\
224 {\\partial f_2 \\over \\partial x_1}\n\
225  & {\\partial f_2 \\over \\partial x_2}\n\
226  & \\cdots\n\
227  & {\\partial f_2 \\over \\partial x_N} \\cr\n\
228  \\vdots & \\vdots & \\ddots & \\vdots \\cr\n\
229 {\\partial f_3 \\over \\partial x_1}\n\
230  & {\\partial f_3 \\over \\partial x_2}\n\
231  & \\cdots\n\
232  & {\\partial f_3 \\over \\partial x_N} \\cr}\\right]$$\n\
233 @end tex\n\
234 @ifnottex\n\
235 \n\
236 @example\n\
237 @group\n\
238  | df_1 df_1 df_1 |\n\
239  | ---- ---- ... ---- |\n\
240  | dx_1 dx_2 dx_N |\n\
241  | |\n\
242  | df_2 df_2 df_2 |\n\
243  | ---- ---- ... ---- |\n\
244  df_i | dx_1 dx_2 dx_N |\n\
245 jac = ---- = | |\n\
246  dx_j | . . . . |\n\
247  | . . . . |\n\
248  | . . . . |\n\
249  | |\n\
250  | df_N df_N df_N |\n\
251  | ---- ---- ... ---- |\n\
252  | dx_1 dx_2 dx_N |\n\
253 @end group\n\
254 @end example\n\
255 \n\
256 @end ifnottex\n\
257 \n\
258 The second and third arguments specify the initial state of the system,\n\
259 @math{x_0}, and the initial value of the independent variable @math{t_0}.\n\
260 \n\
261 The fourth argument is optional, and may be used to specify a set of\n\
262 times that the ODE solver should not integrate past. It is useful for\n\
263 avoiding difficulties with singularities and points where there is a\n\
264 discontinuity in the derivative.\n\
265 \n\
266 After a successful computation, the value of @var{istate} will be 2\n\
267 (consistent with the Fortran version of @sc{lsode}).\n\
268 \n\
269 If the computation is not successful, @var{istate} will be something\n\
270 other than 2 and @var{msg} will contain additional information.\n\
271 \n\
272 You can use the function @code{lsode_options} to set optional\n\
273 parameters for @code{lsode}.\n\
274 @seealso{daspk, dassl, dasrt}\n\
275 @end deftypefn")
276 {
277  octave_value_list retval;
278 
279  warned_fcn_imaginary = false;
280  warned_jac_imaginary = false;
281 
282  unwind_protect frame;
283 
284  frame.protect_var (call_depth);
285  call_depth++;
286 
287  if (call_depth > 1)
288  LSODE_ABORT1 ("invalid recursive call");
289 
290  int nargin = args.length ();
291 
292  if (nargin > 2 && nargin < 5 && nargout < 4)
293  {
294  std::string fcn_name, fname, jac_name, jname;
295  lsode_fcn = 0;
296  lsode_jac = 0;
297 
298  octave_value f_arg = args(0);
299 
300  if (f_arg.is_cell ())
301  {
302  Cell c = f_arg.cell_value ();
303  if (c.length () == 1)
304  f_arg = c(0);
305  else if (c.length () == 2)
306  {
307  if (c(0).is_function_handle () || c(0).is_inline_function ())
308  lsode_fcn = c(0).function_value ();
309  else
310  {
311  fcn_name = unique_symbol_name ("__lsode_fcn__");
312  fname = "function y = ";
313  fname.append (fcn_name);
314  fname.append (" (x, t) y = ");
315  lsode_fcn = extract_function (c(0), "lsode", fcn_name, fname,
316  "; endfunction");
317  }
318 
319  if (lsode_fcn)
320  {
321  if (c(1).is_function_handle () || c(1).is_inline_function ())
322  lsode_jac = c(1).function_value ();
323  else
324  {
325  jac_name = unique_symbol_name ("__lsode_jac__");
326  jname = "function jac = ";
327  jname.append (jac_name);
328  jname.append (" (x, t) jac = ");
329  lsode_jac = extract_function (c(1), "lsode", jac_name,
330  jname, "; endfunction");
331 
332  if (!lsode_jac)
333  {
334  if (fcn_name.length ())
335  clear_function (fcn_name);
336  lsode_fcn = 0;
337  }
338  }
339  }
340  }
341  else
342  LSODE_ABORT1 ("incorrect number of elements in cell array");
343  }
344 
345  if (!lsode_fcn && ! f_arg.is_cell ())
346  {
347  if (f_arg.is_function_handle () || f_arg.is_inline_function ())
348  lsode_fcn = f_arg.function_value ();
349  else
350  {
351  switch (f_arg.rows ())
352  {
353  case 1:
354  do
355  {
356  fcn_name = unique_symbol_name ("__lsode_fcn__");
357  fname = "function y = ";
358  fname.append (fcn_name);
359  fname.append (" (x, t) y = ");
360  lsode_fcn = extract_function (f_arg, "lsode", fcn_name,
361  fname, "; endfunction");
362  }
363  while (0);
364  break;
365 
366  case 2:
367  {
368  string_vector tmp = f_arg.all_strings ();
369 
370  if (! error_state)
371  {
372  fcn_name = unique_symbol_name ("__lsode_fcn__");
373  fname = "function y = ";
374  fname.append (fcn_name);
375  fname.append (" (x, t) y = ");
376  lsode_fcn = extract_function (tmp(0), "lsode", fcn_name,
377  fname, "; endfunction");
378 
379  if (lsode_fcn)
380  {
381  jac_name = unique_symbol_name ("__lsode_jac__");
382  jname = "function jac = ";
383  jname.append (jac_name);
384  jname.append (" (x, t) jac = ");
385  lsode_jac = extract_function (tmp(1), "lsode",
386  jac_name, jname,
387  "; endfunction");
388 
389  if (!lsode_jac)
390  {
391  if (fcn_name.length ())
392  clear_function (fcn_name);
393  lsode_fcn = 0;
394  }
395  }
396  }
397  }
398  break;
399 
400  default:
402  ("first arg should be a string or 2-element string array");
403  }
404  }
405  }
406 
407  if (error_state || ! lsode_fcn)
408  LSODE_ABORT ();
409 
410  ColumnVector state (args(1).vector_value ());
411 
412  if (error_state)
413  LSODE_ABORT1 ("expecting state vector as second argument");
414 
415  ColumnVector out_times (args(2).vector_value ());
416 
417  if (error_state)
418  LSODE_ABORT1 ("expecting output time vector as third argument");
419 
420  ColumnVector crit_times;
421 
422  int crit_times_set = 0;
423  if (nargin > 3)
424  {
425  crit_times = ColumnVector (args(3).vector_value ());
426 
427  if (error_state)
428  LSODE_ABORT1 ("expecting critical time vector as fourth argument");
429 
430  crit_times_set = 1;
431  }
432 
433  double tzero = out_times (0);
434 
436  if (lsode_jac)
438 
439  LSODE ode (state, tzero, func);
440 
441  ode.set_options (lsode_opts);
442 
443  Matrix output;
444  if (crit_times_set)
445  output = ode.integrate (out_times, crit_times);
446  else
447  output = ode.integrate (out_times);
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 = ode.error_message ();
457 
458  retval(2) = msg;
459  retval(1) = static_cast<double> (ode.integration_state ());
460 
461  if (ode.integration_ok ())
462  retval(0) = output;
463  else
464  {
465  retval(0) = Matrix ();
466 
467  if (nargout < 2)
468  error ("lsode: %s", msg.c_str ());
469  }
470  }
471  }
472  else
473  print_usage ();
474 
475  return retval;
476 }
477 
478 /*
479 
480 ## dassl-1.m
481 ##
482 ## Test lsode() function
483 ##
484 ## Author: David Billinghurst ([email protected])
485 ## Comalco Research and Technology
486 ## 20 May 1998
487 ##
488 ## Problem
489 ##
490 ## y1' = -y2, y1(0) = 1
491 ## y2' = y1, y2(0) = 0
492 ##
493 ## Solution
494 ##
495 ## y1(t) = cos(t)
496 ## y2(t) = sin(t)
497 ##
498 %!function xdot = __f (x, t)
499 %! xdot = [-x(2); x(1)];
500 %!endfunction
501 %!test
502 %!
503 %! x0 = [1; 0];
504 %! xdot0 = [0; 1];
505 %! t = (0:1:10)';
506 %!
507 %! tol = 500 * lsode_options ("relative tolerance");
508 %!
509 %! x = lsode ("__f", x0, t);
510 %!
511 %! y = [cos(t), sin(t)];
512 %!
513 %! assert (x, y, tol);
514 
515 %!function xdotdot = __f (x, t)
516 %! xdotdot = [x(2); -x(1)];
517 %!endfunction
518 %!test
519 %!
520 %! x0 = [1; 0];
521 %! t = [0; 2*pi];
522 %! tol = 100 * dassl_options ("relative tolerance");
523 %!
524 %! x = lsode ("__f", x0, t);
525 %!
526 %! y = [1, 0; 1, 0];
527 %!
528 %! assert (x, y, tol);
529 
530 %!function xdot = __f (x, t)
531 %! xdot = x;
532 %!endfunction
533 %!test
534 %!
535 %! x0 = 1;
536 %! t = [0; 1];
537 %! tol = 100 * dassl_options ("relative tolerance");
538 %!
539 %! x = lsode ("__f", x0, t);
540 %!
541 %! y = [1; e];
542 %!
543 %! assert (x, y, tol);
544 
545 %!test
546 %! lsode_options ("absolute tolerance", eps);
547 %! assert (lsode_options ("absolute tolerance") == eps);
548 
549 %!error lsode_options ("foo", 1, 2)
550 */
static int call_depth
Definition: lsode.cc:60
static bool warned_jac_imaginary
Definition: lsode.cc:57
Definition: Cell.h:35
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
Matrix lsode_user_jacobian(const ColumnVector &x, double t)
Definition: lsode.cc:102
virtual ColumnVector integrate(double tt)
Definition: ODE.h:73
void protect_var(T &var)
static bool warned_fcn_imaginary
Definition: lsode.cc:56
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
void error(const char *fmt,...)
Definition: error.cc:476
bool is_cell(void) const
Definition: ov.h:529
ColumnVector lsode_user_function(const ColumnVector &x, double t)
Definition: lsode.cc:63
Definition: LSODE.h:31
void set_options(const LSODE_options &opt)
Definition: LSODE-opts.h:78
bool is_function_handle(void) const
Definition: ov.h:686
string_vector all_strings(bool pad=false) const
Definition: ov.h:894
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
ODEFunc & set_jacobian_function(ODEJacFunc j)
Definition: ODEFunc.h:73
#define LSODE_ABORT()
Definition: lsode.cc:140
static octave_function * lsode_jac
Definition: lsode.cc:53
int error_state
Definition: error.cc:101
bool is_inline_function(void) const
Definition: ov.h:692
Definition: dMatrix.h:35
#define LSODE_ABORT1(msg)
Definition: lsode.cc:143
octave_function * function_value(bool silent=false) const
Definition: ov.cc:1597
void warning(const char *fmt,...)
Definition: error.cc:681
static octave_function * lsode_fcn
Definition: lsode.cc:50
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
static LSODE_options lsode_opts
Definition: LSODE-opts.cc:20
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
void gripe_user_supplied_eval(const char *name)
Definition: gripes.cc:87
F77_RET_T const double * x
std::string error_message(void) const
Definition: LSODE.cc:316