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
lsode.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 "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 Solve the set of differential equations\n\
164 @tex\n\
165 $$ {dx \\over dt} = f (x, t) $$\n\
166 with\n\
167 $$ x(t_0) = x_0 $$\n\
168 @end tex\n\
169 @ifnottex\n\
170 \n\
171 @example\n\
172 @group\n\
173 dx\n\
174 -- = f (x, t)\n\
175 dt\n\
176 @end group\n\
177 @end example\n\
178 \n\
179 @noindent\n\
180 with\n\
181 \n\
182 @example\n\
183 x(t_0) = x_0\n\
184 @end example\n\
185 \n\
186 @end ifnottex\n\
187 The solution is returned in the matrix @var{x}, with each row\n\
188 corresponding to an element of the vector @var{t}. The first element\n\
189 of @var{t} should be @math{t_0} and should correspond to the initial\n\
190 state of the system @var{x_0}, so that the first row of the output\n\
191 is @var{x_0}.\n\
192 \n\
193 The first argument, @var{fcn}, is a string, inline, or function handle\n\
194 that names the function @math{f} to call to compute the vector of right\n\
195 hand sides for the set of equations. The function must have the form\n\
196 \n\
197 @example\n\
198 @var{xdot} = f (@var{x}, @var{t})\n\
199 @end example\n\
200 \n\
201 @noindent\n\
202 in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar.\n\
203 \n\
204 If @var{fcn} is a two-element string array or a two-element cell array\n\
205 of strings, inline functions, or function handles, the first element names\n\
206 the function @math{f} described above, and the second element names a\n\
207 function to compute the Jacobian of @math{f}. The Jacobian function\n\
208 must have the form\n\
209 \n\
210 @example\n\
211 @var{jac} = j (@var{x}, @var{t})\n\
212 @end example\n\
213 \n\
214 @noindent\n\
215 in which @var{jac} is the matrix of partial derivatives\n\
216 @tex\n\
217 $$ J = {\\partial f_i \\over \\partial x_j} = \\left[\\matrix{\n\
218 {\\partial f_1 \\over \\partial x_1}\n\
219  & {\\partial f_1 \\over \\partial x_2}\n\
220  & \\cdots\n\
221  & {\\partial f_1 \\over \\partial x_N} \\cr\n\
222 {\\partial f_2 \\over \\partial x_1}\n\
223  & {\\partial f_2 \\over \\partial x_2}\n\
224  & \\cdots\n\
225  & {\\partial f_2 \\over \\partial x_N} \\cr\n\
226  \\vdots & \\vdots & \\ddots & \\vdots \\cr\n\
227 {\\partial f_3 \\over \\partial x_1}\n\
228  & {\\partial f_3 \\over \\partial x_2}\n\
229  & \\cdots\n\
230  & {\\partial f_3 \\over \\partial x_N} \\cr}\\right]$$\n\
231 @end tex\n\
232 @ifnottex\n\
233 \n\
234 @example\n\
235 @group\n\
236  | df_1 df_1 df_1 |\n\
237  | ---- ---- ... ---- |\n\
238  | dx_1 dx_2 dx_N |\n\
239  | |\n\
240  | df_2 df_2 df_2 |\n\
241  | ---- ---- ... ---- |\n\
242  df_i | dx_1 dx_2 dx_N |\n\
243 jac = ---- = | |\n\
244  dx_j | . . . . |\n\
245  | . . . . |\n\
246  | . . . . |\n\
247  | |\n\
248  | df_N df_N df_N |\n\
249  | ---- ---- ... ---- |\n\
250  | dx_1 dx_2 dx_N |\n\
251 @end group\n\
252 @end example\n\
253 \n\
254 @end ifnottex\n\
255 \n\
256 The second and third arguments specify the initial state of the system,\n\
257 @math{x_0}, and the initial value of the independent variable @math{t_0}.\n\
258 \n\
259 The fourth argument is optional, and may be used to specify a set of\n\
260 times that the ODE 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 2\n\
265 (consistent with the Fortran version of @sc{lsode}).\n\
266 \n\
267 If the computation is not successful, @var{istate} will be something\n\
268 other than 2 and @var{msg} will contain additional information.\n\
269 \n\
270 You can use the function @code{lsode_options} to set optional\n\
271 parameters for @code{lsode}.\n\
272 @seealso{daspk, dassl, dasrt}\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  LSODE_ABORT1 ("invalid recursive call");
287 
288  int nargin = args.length ();
289 
290  if (nargin > 2 && nargin < 5 && nargout < 4)
291  {
292  std::string fcn_name, fname, jac_name, jname;
293  lsode_fcn = 0;
294  lsode_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  lsode_fcn = c(0).function_value ();
307  else
308  {
309  fcn_name = unique_symbol_name ("__lsode_fcn__");
310  fname = "function y = ";
311  fname.append (fcn_name);
312  fname.append (" (x, t) y = ");
313  lsode_fcn = extract_function (c(0), "lsode", fcn_name, fname,
314  "; endfunction");
315  }
316 
317  if (lsode_fcn)
318  {
319  if (c(1).is_function_handle () || c(1).is_inline_function ())
320  lsode_jac = c(1).function_value ();
321  else
322  {
323  jac_name = unique_symbol_name ("__lsode_jac__");
324  jname = "function jac = ";
325  jname.append (jac_name);
326  jname.append (" (x, t) jac = ");
327  lsode_jac = extract_function (c(1), "lsode", jac_name,
328  jname, "; endfunction");
329 
330  if (!lsode_jac)
331  {
332  if (fcn_name.length ())
333  clear_function (fcn_name);
334  lsode_fcn = 0;
335  }
336  }
337  }
338  }
339  else
340  LSODE_ABORT1 ("incorrect number of elements in cell array");
341  }
342 
343  if (!lsode_fcn && ! f_arg.is_cell ())
344  {
345  if (f_arg.is_function_handle () || f_arg.is_inline_function ())
346  lsode_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 ("__lsode_fcn__");
355  fname = "function y = ";
356  fname.append (fcn_name);
357  fname.append (" (x, t) y = ");
358  lsode_fcn = extract_function (f_arg, "lsode", 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 ("__lsode_fcn__");
371  fname = "function y = ";
372  fname.append (fcn_name);
373  fname.append (" (x, t) y = ");
374  lsode_fcn = extract_function (tmp(0), "lsode", fcn_name,
375  fname, "; endfunction");
376 
377  if (lsode_fcn)
378  {
379  jac_name = unique_symbol_name ("__lsode_jac__");
380  jname = "function jac = ";
381  jname.append (jac_name);
382  jname.append (" (x, t) jac = ");
383  lsode_jac = extract_function (tmp(1), "lsode",
384  jac_name, jname,
385  "; endfunction");
386 
387  if (!lsode_jac)
388  {
389  if (fcn_name.length ())
390  clear_function (fcn_name);
391  lsode_fcn = 0;
392  }
393  }
394  }
395  }
396  break;
397 
398  default:
400  ("first arg should be a string or 2-element string array");
401  }
402  }
403  }
404 
405  if (error_state || ! lsode_fcn)
406  LSODE_ABORT ();
407 
408  ColumnVector state (args(1).vector_value ());
409 
410  if (error_state)
411  LSODE_ABORT1 ("expecting state vector as second argument");
412 
413  ColumnVector out_times (args(2).vector_value ());
414 
415  if (error_state)
416  LSODE_ABORT1 ("expecting output time vector as third argument");
417 
418  ColumnVector crit_times;
419 
420  int crit_times_set = 0;
421  if (nargin > 3)
422  {
423  crit_times = ColumnVector (args(3).vector_value ());
424 
425  if (error_state)
426  LSODE_ABORT1 ("expecting critical time vector as fourth argument");
427 
428  crit_times_set = 1;
429  }
430 
431  double tzero = out_times (0);
432 
434  if (lsode_jac)
436 
437  LSODE ode (state, tzero, func);
438 
439  ode.set_options (lsode_opts);
440 
441  Matrix output;
442  if (crit_times_set)
443  output = ode.integrate (out_times, crit_times);
444  else
445  output = ode.integrate (out_times);
446 
447  if (fcn_name.length ())
448  clear_function (fcn_name);
449  if (jac_name.length ())
450  clear_function (jac_name);
451 
452  if (! error_state)
453  {
454  std::string msg = ode.error_message ();
455 
456  retval(2) = msg;
457  retval(1) = static_cast<double> (ode.integration_state ());
458 
459  if (ode.integration_ok ())
460  retval(0) = output;
461  else
462  {
463  retval(0) = Matrix ();
464 
465  if (nargout < 2)
466  error ("lsode: %s", msg.c_str ());
467  }
468  }
469  }
470  else
471  print_usage ();
472 
473  return retval;
474 }
475 
476 /*
477 
478 ## dassl-1.m
479 ##
480 ## Test lsode() function
481 ##
482 ## Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
483 ## Comalco Research and Technology
484 ## 20 May 1998
485 ##
486 ## Problem
487 ##
488 ## y1' = -y2, y1(0) = 1
489 ## y2' = y1, y2(0) = 0
490 ##
491 ## Solution
492 ##
493 ## y1(t) = cos(t)
494 ## y2(t) = sin(t)
495 ##
496 %!function xdot = __f (x, t)
497 %! xdot = [-x(2); x(1)];
498 %!endfunction
499 %!test
500 %!
501 %! x0 = [1; 0];
502 %! xdot0 = [0; 1];
503 %! t = (0:1:10)';
504 %!
505 %! tol = 500 * lsode_options ("relative tolerance");
506 %!
507 %! x = lsode ("__f", x0, t);
508 %!
509 %! y = [cos(t), sin(t)];
510 %!
511 %! assert (x, y, tol);
512 
513 %!function xdotdot = __f (x, t)
514 %! xdotdot = [x(2); -x(1)];
515 %!endfunction
516 %!test
517 %!
518 %! x0 = [1; 0];
519 %! t = [0; 2*pi];
520 %! tol = 100 * dassl_options ("relative tolerance");
521 %!
522 %! x = lsode ("__f", x0, t);
523 %!
524 %! y = [1, 0; 1, 0];
525 %!
526 %! assert (x, y, tol);
527 
528 %!function xdot = __f (x, t)
529 %! xdot = x;
530 %!endfunction
531 %!test
532 %!
533 %! x0 = 1;
534 %! t = [0; 1];
535 %! tol = 100 * dassl_options ("relative tolerance");
536 %!
537 %! x = lsode ("__f", x0, t);
538 %!
539 %! y = [1; e];
540 %!
541 %! assert (x, y, tol);
542 
543 %!test
544 %! lsode_options ("absolute tolerance", eps);
545 %! assert (lsode_options ("absolute tolerance") == eps);
546 
547 %!error lsode_options ("foo", 1, 2)
548 */