lsode.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 1996-2012 John W. Eaton
00004 
00005 This file is part of Octave.
00006 
00007 Octave is free software; you can redistribute it and/or modify it
00008 under the terms of the GNU General Public License as published by the
00009 Free Software Foundation; either version 3 of the License, or (at your
00010 option) any later version.
00011 
00012 Octave is distributed in the hope that it will be useful, but WITHOUT
00013 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00014 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00015 for more details.
00016 
00017 You should have received a copy of the GNU General Public License
00018 along with Octave; see the file COPYING.  If not, see
00019 <http://www.gnu.org/licenses/>.
00020 
00021 */
00022 
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026 
00027 #include <string>
00028 
00029 #include <iomanip>
00030 #include <iostream>
00031 
00032 #include "LSODE.h"
00033 #include "lo-mappers.h"
00034 
00035 #include "defun-dld.h"
00036 #include "error.h"
00037 #include "gripes.h"
00038 #include "oct-obj.h"
00039 #include "ov-fcn.h"
00040 #include "ov-cell.h"
00041 #include "pager.h"
00042 #include "pr-output.h"
00043 #include "unwind-prot.h"
00044 #include "utils.h"
00045 #include "variables.h"
00046 
00047 #include "LSODE-opts.cc"
00048 
00049 // Global pointer for user defined function required by lsode.
00050 static octave_function *lsode_fcn;
00051 
00052 // Global pointer for optional user defined jacobian function used by lsode.
00053 static octave_function *lsode_jac;
00054 
00055 // Have we warned about imaginary values returned from user function?
00056 static bool warned_fcn_imaginary = false;
00057 static bool warned_jac_imaginary = false;
00058 
00059 // Is this a recursive call?
00060 static int call_depth = 0;
00061 
00062 ColumnVector
00063 lsode_user_function (const ColumnVector& x, double t)
00064 {
00065   ColumnVector retval;
00066 
00067   octave_value_list args;
00068   args(1) = t;
00069   args(0) = x;
00070 
00071   if (lsode_fcn)
00072     {
00073       octave_value_list tmp = lsode_fcn->do_multi_index_op (1, args);
00074 
00075       if (error_state)
00076         {
00077           gripe_user_supplied_eval ("lsode");
00078           return retval;
00079         }
00080 
00081       if (tmp.length () > 0 && tmp(0).is_defined ())
00082         {
00083           if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
00084             {
00085               warning ("lsode: ignoring imaginary part returned from user-supplied function");
00086               warned_fcn_imaginary = true;
00087             }
00088 
00089           retval = ColumnVector (tmp(0).vector_value ());
00090 
00091           if (error_state || retval.length () == 0)
00092             gripe_user_supplied_eval ("lsode");
00093         }
00094       else
00095         gripe_user_supplied_eval ("lsode");
00096     }
00097 
00098   return retval;
00099 }
00100 
00101 Matrix
00102 lsode_user_jacobian (const ColumnVector& x, double t)
00103 {
00104   Matrix retval;
00105 
00106   octave_value_list args;
00107   args(1) = t;
00108   args(0) = x;
00109 
00110   if (lsode_jac)
00111     {
00112       octave_value_list tmp = lsode_jac->do_multi_index_op (1, args);
00113 
00114       if (error_state)
00115         {
00116           gripe_user_supplied_eval ("lsode");
00117           return retval;
00118         }
00119 
00120       if (tmp.length () > 0 && tmp(0).is_defined ())
00121         {
00122           if (! warned_jac_imaginary && tmp(0).is_complex_type ())
00123             {
00124               warning ("lsode: ignoring imaginary part returned from user-supplied jacobian function");
00125               warned_jac_imaginary = true;
00126             }
00127 
00128           retval = tmp(0).matrix_value ();
00129 
00130           if (error_state || retval.length () == 0)
00131             gripe_user_supplied_eval ("lsode");
00132         }
00133       else
00134         gripe_user_supplied_eval ("lsode");
00135     }
00136 
00137   return retval;
00138 }
00139 
00140 #define LSODE_ABORT() \
00141   return retval
00142 
00143 #define LSODE_ABORT1(msg) \
00144   do \
00145     { \
00146       ::error ("lsode: " msg); \
00147       LSODE_ABORT (); \
00148     } \
00149   while (0)
00150 
00151 #define LSODE_ABORT2(fmt, arg) \
00152   do \
00153     { \
00154       ::error ("lsode: " fmt, arg); \
00155       LSODE_ABORT (); \
00156     } \
00157   while (0)
00158 
00159 DEFUN_DLD (lsode, args, nargout,
00160   "-*- texinfo -*-\n\
00161 @deftypefn  {Loadable Function} {[@var{x}, @var{istate}, @var{msg}] =} lsode (@var{fcn}, @var{x_0}, @var{t})\n\
00162 @deftypefnx {Loadable Function} {[@var{x}, @var{istate}, @var{msg}] =} lsode (@var{fcn}, @var{x_0}, @var{t}, @var{t_crit})\n\
00163 Solve the set of differential equations\n\
00164 @tex\n\
00165 $$ {dx \\over dt} = f (x, t) $$\n\
00166 with\n\
00167 $$ x(t_0) = x_0 $$\n\
00168 @end tex\n\
00169 @ifnottex\n\
00170 \n\
00171 @example\n\
00172 @group\n\
00173 dx\n\
00174 -- = f(x, t)\n\
00175 dt\n\
00176 @end group\n\
00177 @end example\n\
00178 \n\
00179 @noindent\n\
00180 with\n\
00181 \n\
00182 @example\n\
00183 x(t_0) = x_0\n\
00184 @end example\n\
00185 \n\
00186 @end ifnottex\n\
00187 The solution is returned in the matrix @var{x}, with each row\n\
00188 corresponding to an element of the vector @var{t}.  The first element\n\
00189 of @var{t} should be @math{t_0} and should correspond to the initial\n\
00190 state of the system @var{x_0}, so that the first row of the output\n\
00191 is @var{x_0}.\n\
00192 \n\
00193 The first argument, @var{fcn}, is a string, inline, or function handle\n\
00194 that names the function @math{f} to call to compute the vector of right\n\
00195 hand sides for the set of equations.  The function must have the form\n\
00196 \n\
00197 @example\n\
00198 @var{xdot} = f (@var{x}, @var{t})\n\
00199 @end example\n\
00200 \n\
00201 @noindent\n\
00202 in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar.\n\
00203 \n\
00204 If @var{fcn} is a two-element string array or a two-element cell array\n\
00205 of strings, inline functions, or function handles, the first element names\n\
00206 the function @math{f} described above, and the second element names a\n\
00207 function to compute the Jacobian of @math{f}.  The Jacobian function\n\
00208 must have the form\n\
00209 \n\
00210 @example\n\
00211 @var{jac} = j (@var{x}, @var{t})\n\
00212 @end example\n\
00213 \n\
00214 @noindent\n\
00215 in which @var{jac} is the matrix of partial derivatives\n\
00216 @tex\n\
00217 $$ J = {\\partial f_i \\over \\partial x_j} = \\left[\\matrix{\n\
00218 {\\partial f_1 \\over \\partial x_1}\n\
00219   & {\\partial f_1 \\over \\partial x_2}\n\
00220   & \\cdots\n\
00221   & {\\partial f_1 \\over \\partial x_N} \\cr\n\
00222 {\\partial f_2 \\over \\partial x_1}\n\
00223   & {\\partial f_2 \\over \\partial x_2}\n\
00224   & \\cdots\n\
00225   & {\\partial f_2 \\over \\partial x_N} \\cr\n\
00226  \\vdots & \\vdots & \\ddots & \\vdots \\cr\n\
00227 {\\partial f_3 \\over \\partial x_1}\n\
00228   & {\\partial f_3 \\over \\partial x_2}\n\
00229   & \\cdots\n\
00230   & {\\partial f_3 \\over \\partial x_N} \\cr}\\right]$$\n\
00231 @end tex\n\
00232 @ifnottex\n\
00233 \n\
00234 @example\n\
00235 @group\n\
00236              | df_1  df_1       df_1 |\n\
00237              | ----  ----  ...  ---- |\n\
00238              | dx_1  dx_2       dx_N |\n\
00239              |                       |\n\
00240              | df_2  df_2       df_2 |\n\
00241              | ----  ----  ...  ---- |\n\
00242       df_i   | dx_1  dx_2       dx_N |\n\
00243 jac = ---- = |                       |\n\
00244       dx_j   |  .    .     .    .    |\n\
00245              |  .    .      .   .    |\n\
00246              |  .    .       .  .    |\n\
00247              |                       |\n\
00248              | df_N  df_N       df_N |\n\
00249              | ----  ----  ...  ---- |\n\
00250              | dx_1  dx_2       dx_N |\n\
00251 @end group\n\
00252 @end example\n\
00253 \n\
00254 @end ifnottex\n\
00255 \n\
00256 The second and third arguments specify the initial state of the system,\n\
00257 @math{x_0}, and the initial value of the independent variable @math{t_0}.\n\
00258 \n\
00259 The fourth argument is optional, and may be used to specify a set of\n\
00260 times that the ODE solver should not integrate past.  It is useful for\n\
00261 avoiding difficulties with singularities and points where there is a\n\
00262 discontinuity in the derivative.\n\
00263 \n\
00264 After a successful computation, the value of @var{istate} will be 2\n\
00265 (consistent with the Fortran version of @sc{lsode}).\n\
00266 \n\
00267 If the computation is not successful, @var{istate} will be something\n\
00268 other than 2 and @var{msg} will contain additional information.\n\
00269 \n\
00270 You can use the function @code{lsode_options} to set optional\n\
00271 parameters for @code{lsode}.\n\
00272 @seealso{daspk, dassl, dasrt}\n\
00273 @end deftypefn")
00274 {
00275   octave_value_list retval;
00276 
00277   warned_fcn_imaginary = false;
00278   warned_jac_imaginary = false;
00279 
00280   unwind_protect frame;
00281 
00282   frame.protect_var (call_depth);
00283   call_depth++;
00284 
00285   if (call_depth > 1)
00286     LSODE_ABORT1 ("invalid recursive call");
00287 
00288   int nargin = args.length ();
00289 
00290   if (nargin > 2 && nargin < 5 && nargout < 4)
00291     {
00292       std::string fcn_name, fname, jac_name, jname;
00293       lsode_fcn = 0;
00294       lsode_jac = 0;
00295 
00296       octave_value f_arg = args(0);
00297 
00298       if (f_arg.is_cell ())
00299         {
00300           Cell c = f_arg.cell_value ();
00301           if (c.length() == 1)
00302             f_arg = c(0);
00303           else if (c.length() == 2)
00304             {
00305               if (c(0).is_function_handle () || c(0).is_inline_function ())
00306                 lsode_fcn = c(0).function_value ();
00307               else
00308                 {
00309                   fcn_name = unique_symbol_name ("__lsode_fcn__");
00310                   fname = "function y = ";
00311                   fname.append (fcn_name);
00312                   fname.append (" (x, t) y = ");
00313                   lsode_fcn = extract_function
00314                     (c(0), "lsode", fcn_name, fname, "; endfunction");
00315                 }
00316 
00317               if (lsode_fcn)
00318                 {
00319                   if (c(1).is_function_handle () || c(1).is_inline_function ())
00320                     lsode_jac = c(1).function_value ();
00321                   else
00322                     {
00323                         jac_name = unique_symbol_name ("__lsode_jac__");
00324                         jname = "function jac = ";
00325                         jname.append(jac_name);
00326                         jname.append (" (x, t) jac = ");
00327                         lsode_jac = extract_function
00328                           (c(1), "lsode", jac_name, jname, "; endfunction");
00329 
00330                       if (!lsode_jac)
00331                         {
00332                           if (fcn_name.length())
00333                             clear_function (fcn_name);
00334                           lsode_fcn = 0;
00335                         }
00336                     }
00337                 }
00338             }
00339           else
00340             LSODE_ABORT1 ("incorrect number of elements in cell array");
00341         }
00342 
00343       if (!lsode_fcn && ! f_arg.is_cell())
00344         {
00345           if (f_arg.is_function_handle () || f_arg.is_inline_function ())
00346             lsode_fcn = f_arg.function_value ();
00347           else
00348             {
00349               switch (f_arg.rows ())
00350                 {
00351                 case 1:
00352                   do
00353                     {
00354                       fcn_name = unique_symbol_name ("__lsode_fcn__");
00355                       fname = "function y = ";
00356                       fname.append (fcn_name);
00357                       fname.append (" (x, t) y = ");
00358                       lsode_fcn = extract_function
00359                         (f_arg, "lsode", fcn_name, fname, "; endfunction");
00360                     }
00361                   while (0);
00362                   break;
00363 
00364                 case 2:
00365                   {
00366                     string_vector tmp = f_arg.all_strings ();
00367 
00368                     if (! error_state)
00369                       {
00370                         fcn_name = unique_symbol_name ("__lsode_fcn__");
00371                         fname = "function y = ";
00372                         fname.append (fcn_name);
00373                         fname.append (" (x, t) y = ");
00374                         lsode_fcn = extract_function
00375                           (tmp(0), "lsode", fcn_name, fname, "; endfunction");
00376 
00377                         if (lsode_fcn)
00378                           {
00379                             jac_name = unique_symbol_name ("__lsode_jac__");
00380                             jname = "function jac = ";
00381                             jname.append(jac_name);
00382                             jname.append (" (x, t) jac = ");
00383                             lsode_jac = extract_function
00384                               (tmp(1), "lsode", jac_name, jname,
00385                               "; endfunction");
00386 
00387                             if (!lsode_jac)
00388                               {
00389                                 if (fcn_name.length())
00390                                   clear_function (fcn_name);
00391                                 lsode_fcn = 0;
00392                               }
00393                           }
00394                       }
00395                   }
00396                   break;
00397 
00398                 default:
00399                   LSODE_ABORT1
00400                     ("first arg should be a string or 2-element string array");
00401                 }
00402             }
00403         }
00404 
00405       if (error_state || ! lsode_fcn)
00406         LSODE_ABORT ();
00407 
00408       ColumnVector state (args(1).vector_value ());
00409 
00410       if (error_state)
00411         LSODE_ABORT1 ("expecting state vector as second argument");
00412 
00413       ColumnVector out_times (args(2).vector_value ());
00414 
00415       if (error_state)
00416         LSODE_ABORT1 ("expecting output time vector as third argument");
00417 
00418       ColumnVector crit_times;
00419 
00420       int crit_times_set = 0;
00421       if (nargin > 3)
00422         {
00423           crit_times = ColumnVector (args(3).vector_value ());
00424 
00425           if (error_state)
00426             LSODE_ABORT1 ("expecting critical time vector as fourth argument");
00427 
00428           crit_times_set = 1;
00429         }
00430 
00431       double tzero = out_times (0);
00432 
00433       ODEFunc func (lsode_user_function);
00434       if (lsode_jac)
00435         func.set_jacobian_function (lsode_user_jacobian);
00436 
00437       LSODE ode (state, tzero, func);
00438 
00439       ode.set_options (lsode_opts);
00440 
00441       Matrix output;
00442       if (crit_times_set)
00443         output = ode.integrate (out_times, crit_times);
00444       else
00445         output = ode.integrate (out_times);
00446 
00447       if (fcn_name.length())
00448         clear_function (fcn_name);
00449       if (jac_name.length())
00450         clear_function (jac_name);
00451 
00452       if (! error_state)
00453         {
00454           std::string msg = ode.error_message ();
00455 
00456           retval(2) = msg;
00457           retval(1) = static_cast<double> (ode.integration_state ());
00458 
00459           if (ode.integration_ok ())
00460             retval(0) = output;
00461           else
00462             {
00463               retval(0) = Matrix ();
00464 
00465               if (nargout < 2)
00466                 error ("lsode: %s", msg.c_str ());
00467             }
00468         }
00469     }
00470   else
00471     print_usage ();
00472 
00473   return retval;
00474 }
00475 
00476 /*
00477 
00478 %% dassl-1.m
00479 %%
00480 %% Test lsode() function
00481 %%
00482 %% Author: David Billinghurst (David.Billinghurst@riotinto.com.au)
00483 %%         Comalco Research and Technology
00484 %%         20 May 1998
00485 %%
00486 %% Problem
00487 %%
00488 %%    y1' = -y2,   y1(0) = 1
00489 %%    y2' =  y1,   y2(0) = 0
00490 %%
00491 %% Solution
00492 %%
00493 %%    y1(t) = cos(t)
00494 %%    y2(t) = sin(t)
00495 %!function xdot = __f (x, t)
00496 %!  xdot = [-x(2); x(1)];
00497 %!test
00498 %!
00499 %! x0 = [1; 0];
00500 %! xdot0 = [0; 1];
00501 %! t = (0:1:10)';
00502 %!
00503 %! tol = 500 * lsode_options ("relative tolerance");
00504 %!
00505 %!
00506 %! x = lsode ("__f", x0, t);
00507 %!
00508 %! y = [cos(t), sin(t)];
00509 %!
00510 %! assert(all (all (abs (x - y) < tol)));
00511 
00512 %!function xdotdot = __f (x, t)
00513 %!  xdotdot = [x(2); -x(1)];
00514 %!test
00515 %!
00516 %! x0 = [1; 0];
00517 %! t = [0; 2*pi];
00518 %! tol = 100 * dassl_options ("relative tolerance");
00519 %!
00520 %! x = lsode ("__f", x0, t);
00521 %!
00522 %! y = [1, 0; 1, 0];
00523 %!
00524 %! assert(all (all (abs (x - y) < tol)));
00525 
00526 %!function xdot = __f (x, t)
00527 %!  xdot = x;
00528 %!test
00529 %!
00530 %! x0 = 1;
00531 %! t = [0; 1];
00532 %! tol = 100 * dassl_options ("relative tolerance");
00533 %!
00534 %! x = lsode ("__f", x0, t);
00535 %!
00536 %! y = [1; e];
00537 %!
00538 %! assert(all (all (abs (x - y) < tol)));
00539 
00540 %!test
00541 %! lsode_options ("absolute tolerance", eps);
00542 %! assert(lsode_options ("absolute tolerance") == eps);
00543 
00544 %!error <Invalid call to lsode_options> lsode_options ("foo", 1, 2);
00545 
00546 */
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines