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) 1993-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 <cfloat>
28 
29 #include <sstream>
30 
31 #include "LSODE.h"
32 #include "f77-fcn.h"
33 #include "lo-error.h"
34 #include "lo-math.h"
35 #include "quit.h"
36 
38  const double&, double*,
39  double*, octave_idx_type&);
40 
42  const double&, double*,
43  const octave_idx_type&,
44  const octave_idx_type&,
45  double*, const octave_idx_type&);
46 
47 extern "C"
48 {
49  F77_RET_T
51  double&, double&, octave_idx_type&, double&,
52  const double*, octave_idx_type&,
53  octave_idx_type&, octave_idx_type&,
54  double*, octave_idx_type&, octave_idx_type*,
55  octave_idx_type&, lsode_jac_ptr,
56  octave_idx_type&);
57 }
58 
62 
63 static octave_idx_type
64 lsode_f (const octave_idx_type& neq, const double& time, double *,
65  double *deriv, octave_idx_type& ierr)
66 {
67  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
68 
69  ColumnVector tmp_deriv;
70 
71  // NOTE: this won't work if LSODE passes copies of the state vector.
72  // In that case we have to create a temporary vector object
73  // and copy.
74 
75  tmp_deriv = (*user_fun) (*tmp_x, time);
76 
77  if (tmp_deriv.length () == 0)
78  ierr = -1;
79  else
80  {
81  for (octave_idx_type i = 0; i < neq; i++)
82  deriv[i] = tmp_deriv.elem (i);
83  }
84 
85  END_INTERRUPT_WITH_EXCEPTIONS;
86 
87  return 0;
88 }
89 
90 static octave_idx_type
91 lsode_j (const octave_idx_type& neq, const double& time, double *,
92  const octave_idx_type&, const octave_idx_type&, double *pd,
93  const octave_idx_type& nrowpd)
94 {
95  BEGIN_INTERRUPT_WITH_EXCEPTIONS;
96 
97  Matrix tmp_jac (neq, neq);
98 
99  // NOTE: this won't work if LSODE passes copies of the state vector.
100  // In that case we have to create a temporary vector object
101  // and copy.
102 
103  tmp_jac = (*user_jac) (*tmp_x, time);
104 
105  for (octave_idx_type j = 0; j < neq; j++)
106  for (octave_idx_type i = 0; i < neq; i++)
107  pd[nrowpd * j + i] = tmp_jac (i, j);
108 
109  END_INTERRUPT_WITH_EXCEPTIONS;
110 
111  return 0;
112 }
113 
115 LSODE::do_integrate (double tout)
116 {
117  ColumnVector retval;
118 
119  static octave_idx_type nn = 0;
120 
122  {
123  integration_error = false;
124 
125  initialized = true;
126 
127  istate = 1;
128 
129  octave_idx_type n = size ();
130 
131  nn = n;
132 
133  octave_idx_type max_maxord = 0;
134 
135  if (integration_method () == "stiff")
136  {
137  max_maxord = 5;
138 
139  if (jac)
140  method_flag = 21;
141  else
142  method_flag = 22;
143 
144  liw = 20 + n;
145  lrw = 22 + n * (9 + n);
146  }
147  else
148  {
149  max_maxord = 12;
150 
151  method_flag = 10;
152 
153  liw = 20;
154  lrw = 22 + 16 * n;
155  }
156 
157  maxord = maximum_order ();
158 
159  iwork.resize (dim_vector (liw, 1));
160 
161  for (octave_idx_type i = 4; i < 9; i++)
162  iwork(i) = 0;
163 
164  rwork.resize (dim_vector (lrw, 1));
165 
166  for (octave_idx_type i = 4; i < 9; i++)
167  rwork(i) = 0;
168 
169  if (maxord >= 0)
170  {
171  if (maxord > 0 && maxord <= max_maxord)
172  {
173  iwork(4) = maxord;
174  iopt = 1;
175  }
176  else
177  {
178  (*current_liboctave_error_handler)
179  ("lsode: invalid value for maximum order");
180  integration_error = true;
181  return retval;
182  }
183  }
184 
185  if (stop_time_set)
186  {
187  itask = 4;
188  rwork(0) = stop_time;
189  iopt = 1;
190  }
191  else
192  {
193  itask = 1;
194  }
195 
196  restart = false;
197 
198  // ODEFunc
199 
200  // NOTE: this won't work if LSODE passes copies of the state vector.
201  // In that case we have to create a temporary vector object
202  // and copy.
203 
204  tmp_x = &x;
205 
206  user_fun = function ();
208 
209  ColumnVector xdot = (*user_fun) (x, t);
210 
211  if (x.length () != xdot.length ())
212  {
213  (*current_liboctave_error_handler)
214  ("lsode: inconsistent sizes for state and derivative vectors");
215 
216  integration_error = true;
217  return retval;
218  }
219 
220  ODEFunc::reset = false;
221 
222  // LSODE_options
223 
226 
227  octave_idx_type abs_tol_len = abs_tol.length ();
228 
229  if (abs_tol_len == 1)
230  itol = 1;
231  else if (abs_tol_len == n)
232  itol = 2;
233  else
234  {
235  (*current_liboctave_error_handler)
236  ("lsode: inconsistent sizes for state and absolute tolerance vectors");
237 
238  integration_error = true;
239  return retval;
240  }
241 
242  double iss = initial_step_size ();
243  if (iss >= 0.0)
244  {
245  rwork(4) = iss;
246  iopt = 1;
247  }
248 
249  double maxss = maximum_step_size ();
250  if (maxss >= 0.0)
251  {
252  rwork(5) = maxss;
253  iopt = 1;
254  }
255 
256  double minss = minimum_step_size ();
257  if (minss >= 0.0)
258  {
259  rwork(6) = minss;
260  iopt = 1;
261  }
262 
263  octave_idx_type sl = step_limit ();
264  if (sl > 0)
265  {
266  iwork(5) = sl;
267  iopt = 1;
268  }
269 
270  LSODE_options::reset = false;
271  }
272 
273  double *px = x.fortran_vec ();
274 
275  double *pabs_tol = abs_tol.fortran_vec ();
276 
277  octave_idx_type *piwork = iwork.fortran_vec ();
278  double *prwork = rwork.fortran_vec ();
279 
280  F77_XFCN (dlsode, DLSODE, (lsode_f, nn, px, t, tout, itol, rel_tol,
281  pabs_tol, itask, istate, iopt, prwork, lrw,
282  piwork, liw, lsode_j, method_flag));
283 
284  switch (istate)
285  {
286  case 1: // prior to initial integration step.
287  case 2: // lsode was successful.
288  retval = x;
289  t = tout;
290  break;
291 
292  case -1: // excess work done on this call (perhaps wrong mf).
293  case -2: // excess accuracy requested (tolerances too small).
294  case -3: // invalid input detected (see printed message).
295  case -4: // repeated error test failures (check all inputs).
296  case -5: // repeated convergence failures (perhaps bad Jacobian
297  // supplied or wrong choice of mf or tolerances).
298  case -6: // error weight became zero during problem. (solution
299  // component i vanished, and atol or atol(i) = 0.)
300  case -13: // return requested in user-supplied function.
301  integration_error = true;
302  break;
303 
304  default:
305  integration_error = true;
306  (*current_liboctave_error_handler)
307  ("unrecognized value of istate (= %d) returned from lsode",
308  istate);
309  break;
310  }
311 
312  return retval;
313 }
314 
315 std::string
317 {
318  std::string retval;
319 
320  std::ostringstream buf;
321  buf << t;
322  std::string t_curr = buf.str ();
323 
324  switch (istate)
325  {
326  case 1:
327  retval = "prior to initial integration step";
328  break;
329 
330  case 2:
331  retval = "successful exit";
332  break;
333 
334  case 3:
335  retval = "prior to continuation call with modified parameters";
336  break;
337 
338  case -1:
339  retval = std::string ("excess work on this call (t = ")
340  + t_curr + "; perhaps wrong integration method)";
341  break;
342 
343  case -2:
344  retval = "excess accuracy requested (tolerances too small)";
345  break;
346 
347  case -3:
348  retval = "invalid input detected (see printed message)";
349  break;
350 
351  case -4:
352  retval = std::string ("repeated error test failures (t = ")
353  + t_curr + "; check all inputs)";
354  break;
355 
356  case -5:
357  retval = std::string ("repeated convergence failures (t = ")
358  + t_curr
359  + "; perhaps bad Jacobian supplied or wrong choice of integration method or tolerances)";
360  break;
361 
362  case -6:
363  retval = std::string ("error weight became zero during problem. (t = ")
364  + t_curr
365  + "; solution component i vanished, and atol or atol(i) == 0)";
366  break;
367 
368  case -13:
369  retval = "return requested in user-supplied function (t = "
370  + t_curr + ")";
371  break;
372 
373  default:
374  retval = "unknown error state";
375  break;
376  }
377 
378  return retval;
379 }
380 
381 Matrix
383 {
384  Matrix retval;
385 
386  octave_idx_type n_out = tout.capacity ();
387  octave_idx_type n = size ();
388 
389  if (n_out > 0 && n > 0)
390  {
391  retval.resize (n_out, n);
392 
393  for (octave_idx_type i = 0; i < n; i++)
394  retval.elem (0, i) = x.elem (i);
395 
396  for (octave_idx_type j = 1; j < n_out; j++)
397  {
398  ColumnVector x_next = do_integrate (tout.elem (j));
399 
400  if (integration_error)
401  return retval;
402 
403  for (octave_idx_type i = 0; i < n; i++)
404  retval.elem (j, i) = x_next.elem (i);
405  }
406  }
407 
408  return retval;
409 }
410 
411 Matrix
412 LSODE::do_integrate (const ColumnVector& tout, const ColumnVector& tcrit)
413 {
414  Matrix retval;
415 
416  octave_idx_type n_out = tout.capacity ();
417  octave_idx_type n = size ();
418 
419  if (n_out > 0 && n > 0)
420  {
421  retval.resize (n_out, n);
422 
423  for (octave_idx_type i = 0; i < n; i++)
424  retval.elem (0, i) = x.elem (i);
425 
426  octave_idx_type n_crit = tcrit.capacity ();
427 
428  if (n_crit > 0)
429  {
430  octave_idx_type i_crit = 0;
431  octave_idx_type i_out = 1;
432  double next_crit = tcrit.elem (0);
433  double next_out;
434  while (i_out < n_out)
435  {
436  bool do_restart = false;
437 
438  next_out = tout.elem (i_out);
439  if (i_crit < n_crit)
440  next_crit = tcrit.elem (i_crit);
441 
442  octave_idx_type save_output;
443  double t_out;
444 
445  if (next_crit == next_out)
446  {
447  set_stop_time (next_crit);
448  t_out = next_out;
449  save_output = 1;
450  i_out++;
451  i_crit++;
452  do_restart = true;
453  }
454  else if (next_crit < next_out)
455  {
456  if (i_crit < n_crit)
457  {
458  set_stop_time (next_crit);
459  t_out = next_crit;
460  save_output = 0;
461  i_crit++;
462  do_restart = true;
463  }
464  else
465  {
466  clear_stop_time ();
467  t_out = next_out;
468  save_output = 1;
469  i_out++;
470  }
471  }
472  else
473  {
474  set_stop_time (next_crit);
475  t_out = next_out;
476  save_output = 1;
477  i_out++;
478  }
479 
480  ColumnVector x_next = do_integrate (t_out);
481 
482  if (integration_error)
483  return retval;
484 
485  if (save_output)
486  {
487  for (octave_idx_type i = 0; i < n; i++)
488  retval.elem (i_out-1, i) = x_next.elem (i);
489  }
490 
491  if (do_restart)
492  force_restart ();
493  }
494  }
495  else
496  {
497  retval = do_integrate (tout);
498 
499  if (integration_error)
500  return retval;
501  }
502  }
503 
504  return retval;
505 }