dasrt.cc

Go to the documentation of this file.
00001 /*
00002 
00003 Copyright (C) 2002-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 <iostream>
00028 #include <string>
00029 
00030 #include "DASRT.h"
00031 #include "lo-mappers.h"
00032 
00033 #include "defun-dld.h"
00034 #include "error.h"
00035 #include "gripes.h"
00036 #include "oct-obj.h"
00037 #include "ov-fcn.h"
00038 #include "ov-cell.h"
00039 #include "pager.h"
00040 #include "parse.h"
00041 #include "unwind-prot.h"
00042 #include "utils.h"
00043 #include "variables.h"
00044 
00045 #include "DASRT-opts.cc"
00046 
00047 // Global pointers for user defined function required by dasrt.
00048 static octave_function *dasrt_f;
00049 static octave_function *dasrt_j;
00050 static octave_function *dasrt_cf;
00051 
00052 // Have we warned about imaginary values returned from user function?
00053 static bool warned_fcn_imaginary = false;
00054 static bool warned_jac_imaginary = false;
00055 static bool warned_cf_imaginary = false;
00056 
00057 // Is this a recursive call?
00058 static int call_depth = 0;
00059 
00060 static ColumnVector
00061 dasrt_user_f (const ColumnVector& x, const ColumnVector& xdot,
00062               double t, octave_idx_type&)
00063 {
00064   ColumnVector retval;
00065 
00066   assert (x.capacity () == xdot.capacity ());
00067 
00068   octave_value_list args;
00069 
00070   args(2) = t;
00071   args(1) = xdot;
00072   args(0) = x;
00073 
00074   if (dasrt_f)
00075     {
00076       octave_value_list tmp = dasrt_f->do_multi_index_op (1, args);
00077 
00078       if (error_state)
00079         {
00080           gripe_user_supplied_eval ("dasrt");
00081           return retval;
00082         }
00083 
00084       if (tmp.length () > 0 && tmp(0).is_defined ())
00085         {
00086           if (! warned_fcn_imaginary && tmp(0).is_complex_type ())
00087             {
00088               warning ("dasrt: ignoring imaginary part returned from user-supplied function");
00089               warned_fcn_imaginary = true;
00090             }
00091 
00092           retval = ColumnVector (tmp(0).vector_value ());
00093 
00094           if (error_state || retval.length () == 0)
00095             gripe_user_supplied_eval ("dasrt");
00096         }
00097       else
00098         gripe_user_supplied_eval ("dasrt");
00099     }
00100 
00101   return retval;
00102 }
00103 
00104 static ColumnVector
00105 dasrt_user_cf (const ColumnVector& x, double t)
00106 {
00107   ColumnVector retval;
00108 
00109   octave_value_list args;
00110 
00111   args(1) = t;
00112   args(0) = x;
00113 
00114   if (dasrt_cf)
00115     {
00116       octave_value_list tmp = dasrt_cf->do_multi_index_op (1, args);
00117 
00118       if (error_state)
00119         {
00120           gripe_user_supplied_eval ("dasrt");
00121           return retval;
00122         }
00123 
00124       if (tmp.length () > 0 && tmp(0).is_defined ())
00125         {
00126           if (! warned_cf_imaginary && tmp(0).is_complex_type ())
00127             {
00128               warning ("dasrt: ignoring imaginary part returned from user-supplied constraint function");
00129               warned_cf_imaginary = true;
00130             }
00131 
00132           retval = ColumnVector (tmp(0).vector_value ());
00133 
00134           if (error_state || retval.length () == 0)
00135             gripe_user_supplied_eval ("dasrt");
00136         }
00137       else
00138         gripe_user_supplied_eval ("dasrt");
00139     }
00140 
00141   return retval;
00142 }
00143 
00144 static Matrix
00145 dasrt_user_j (const ColumnVector& x, const ColumnVector& xdot,
00146               double t, double cj)
00147 {
00148   Matrix retval;
00149 
00150   assert (x.capacity () == xdot.capacity ());
00151 
00152   octave_value_list args;
00153 
00154   args(3) = cj;
00155   args(2) = t;
00156   args(1) = xdot;
00157   args(0) = x;
00158 
00159   if (dasrt_j)
00160     {
00161       octave_value_list tmp = dasrt_j->do_multi_index_op (1, args);
00162 
00163       if (error_state)
00164         {
00165           gripe_user_supplied_eval ("dasrt");
00166           return retval;
00167         }
00168 
00169       int tlen = tmp.length ();
00170       if (tlen > 0 && tmp(0).is_defined ())
00171         {
00172           if (! warned_jac_imaginary && tmp(0).is_complex_type ())
00173             {
00174               warning ("dasrt: ignoring imaginary part returned from user-supplied jacobian function");
00175               warned_jac_imaginary = true;
00176             }
00177 
00178           retval = tmp(0).matrix_value ();
00179 
00180           if (error_state || retval.length () == 0)
00181             gripe_user_supplied_eval ("dasrt");
00182         }
00183       else
00184         gripe_user_supplied_eval ("dasrt");
00185     }
00186 
00187   return retval;
00188 }
00189 
00190 #define DASRT_ABORT \
00191   return retval
00192 
00193 #define DASRT_ABORT1(msg) \
00194   do \
00195     { \
00196       ::error ("dasrt: " msg); \
00197       DASRT_ABORT; \
00198     } \
00199   while (0)
00200 
00201 #define DASRT_ABORT2(fmt, arg) \
00202   do \
00203     { \
00204       ::error ("dasrt: " fmt, arg); \
00205       DASRT_ABORT; \
00206     } \
00207   while (0)
00208 
00209 DEFUN_DLD (dasrt, args, nargout,
00210   "-*- texinfo -*-\n\
00211 @deftypefn  {Loadable Function} {[@var{x}, @var{xdot}, @var{t_out}, @var{istat}, @var{msg}] =} dasrt (@var{fcn}, [], @var{x_0}, @var{xdot_0}, @var{t})\n\
00212 @deftypefnx {Loadable Function} {@dots{} =} dasrt (@var{fcn}, @var{g}, @var{x_0}, @var{xdot_0}, @var{t})\n\
00213 @deftypefnx {Loadable Function} {@dots{} =} dasrt (@var{fcn}, [], @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})\n\
00214 @deftypefnx {Loadable Function} {@dots{} =} dasrt (@var{fcn}, @var{g}, @var{x_0}, @var{xdot_0}, @var{t}, @var{t_crit})\n\
00215 Solve the set of differential-algebraic equations\n\
00216 @tex\n\
00217 $$ 0 = f (x, \\dot{x}, t) $$\n\
00218 with\n\
00219 $$ x(t_0) = x_0, \\dot{x}(t_0) = \\dot{x}_0 $$\n\
00220 @end tex\n\
00221 @ifnottex\n\
00222 \n\
00223 @example\n\
00224 0 = f (x, xdot, t)\n\
00225 @end example\n\
00226 \n\
00227 @noindent\n\
00228 with\n\
00229 \n\
00230 @example\n\
00231 x(t_0) = x_0, xdot(t_0) = xdot_0\n\
00232 @end example\n\
00233 \n\
00234 @end ifnottex\n\
00235 with functional stopping criteria (root solving).\n\
00236 \n\
00237 The solution is returned in the matrices @var{x} and @var{xdot},\n\
00238 with each row in the result matrices corresponding to one of the\n\
00239 elements in the vector @var{t_out}.  The first element of @var{t}\n\
00240 should be @math{t_0} and correspond to the initial state of the\n\
00241 system @var{x_0} and its derivative @var{xdot_0}, so that the first\n\
00242 row of the output @var{x} is @var{x_0} and the first row\n\
00243 of the output @var{xdot} is @var{xdot_0}.\n\
00244 \n\
00245 The vector @var{t} provides an upper limit on the length of the\n\
00246 integration.  If the stopping condition is met, the vector\n\
00247 @var{t_out} will be shorter than @var{t}, and the final element of\n\
00248 @var{t_out} will be the point at which the stopping condition was met,\n\
00249 and may not correspond to any element of the vector @var{t}.\n\
00250 \n\
00251 The first argument, @var{fcn}, is a string, inline, or function handle\n\
00252 that names the function @math{f} to call to compute the vector of\n\
00253 residuals for the set of equations.  It must have the form\n\
00254 \n\
00255 @example\n\
00256 @var{res} = f (@var{x}, @var{xdot}, @var{t})\n\
00257 @end example\n\
00258 \n\
00259 @noindent\n\
00260 in which @var{x}, @var{xdot}, and @var{res} are vectors, and @var{t} is a\n\
00261 scalar.\n\
00262 \n\
00263 If @var{fcn} is a two-element string array or a two-element cell array\n\
00264 of strings, inline functions, or function handles, the first element names\n\
00265 the function @math{f} described above, and the second element names a\n\
00266 function to compute the modified Jacobian\n\
00267 \n\
00268 @tex\n\
00269 $$\n\
00270 J = {\\partial f \\over \\partial x}\n\
00271   + c {\\partial f \\over \\partial \\dot{x}}\n\
00272 $$\n\
00273 @end tex\n\
00274 @ifnottex\n\
00275 \n\
00276 @example\n\
00277 @group\n\
00278       df       df\n\
00279 jac = -- + c ------\n\
00280       dx     d xdot\n\
00281 @end group\n\
00282 @end example\n\
00283 \n\
00284 @end ifnottex\n\
00285 \n\
00286 The modified Jacobian function must have the form\n\
00287 \n\
00288 @example\n\
00289 @group\n\
00290 \n\
00291 @var{jac} = j (@var{x}, @var{xdot}, @var{t}, @var{c})\n\
00292 \n\
00293 @end group\n\
00294 @end example\n\
00295 \n\
00296 The optional second argument names a function that defines the\n\
00297 constraint functions whose roots are desired during the integration.\n\
00298 This function must have the form\n\
00299 \n\
00300 @example\n\
00301 @var{g_out} = g (@var{x}, @var{t})\n\
00302 @end example\n\
00303 \n\
00304 @noindent\n\
00305 and return a vector of the constraint function values.\n\
00306 If the value of any of the constraint functions changes sign, @sc{dasrt}\n\
00307 will attempt to stop the integration at the point of the sign change.\n\
00308 \n\
00309 If the name of the constraint function is omitted, @code{dasrt} solves\n\
00310 the same problem as @code{daspk} or @code{dassl}.\n\
00311 \n\
00312 Note that because of numerical errors in the constraint functions\n\
00313 due to round-off and integration error, @sc{dasrt} may return false\n\
00314 roots, or return the same root at two or more nearly equal values of\n\
00315 @var{T}.  If such false roots are suspected, the user should consider\n\
00316 smaller error tolerances or higher precision in the evaluation of the\n\
00317 constraint functions.\n\
00318 \n\
00319 If a root of some constraint function defines the end of the problem,\n\
00320 the input to @sc{dasrt} should nevertheless allow integration to a\n\
00321 point slightly past that root, so that @sc{dasrt} can locate the root\n\
00322 by interpolation.\n\
00323 \n\
00324 The third and fourth arguments to @code{dasrt} specify the initial\n\
00325 condition of the states and their derivatives, and the fourth argument\n\
00326 specifies a vector of output times at which the solution is desired,\n\
00327 including the time corresponding to the initial condition.\n\
00328 \n\
00329 The set of initial states and derivatives are not strictly required to\n\
00330 be consistent.  In practice, however, @sc{dassl} is not very good at\n\
00331 determining a consistent set for you, so it is best if you ensure that\n\
00332 the initial values result in the function evaluating to zero.\n\
00333 \n\
00334 The sixth argument is optional, and may be used to specify a set of\n\
00335 times that the DAE solver should not integrate past.  It is useful for\n\
00336 avoiding difficulties with singularities and points where there is a\n\
00337 discontinuity in the derivative.\n\
00338 \n\
00339 After a successful computation, the value of @var{istate} will be\n\
00340 greater than zero (consistent with the Fortran version of @sc{dassl}).\n\
00341 \n\
00342 If the computation is not successful, the value of @var{istate} will be\n\
00343 less than zero and @var{msg} will contain additional information.\n\
00344 \n\
00345 You can use the function @code{dasrt_options} to set optional\n\
00346 parameters for @code{dasrt}.\n\
00347 @seealso{dasrt_options, daspk, dasrt, lsode}\n\
00348 @end deftypefn")
00349 {
00350   octave_value_list retval;
00351 
00352   warned_fcn_imaginary = false;
00353   warned_jac_imaginary = false;
00354   warned_cf_imaginary = false;
00355 
00356   unwind_protect frame;
00357 
00358   frame.protect_var (call_depth);
00359   call_depth++;
00360 
00361   if (call_depth > 1)
00362     DASRT_ABORT1 ("invalid recursive call");
00363 
00364   int argp = 0;
00365 
00366   int nargin = args.length ();
00367 
00368   if (nargin < 4 || nargin > 6)
00369     {
00370       print_usage ();
00371       return retval;
00372     }
00373 
00374   std::string fcn_name, fname, jac_name, jname;
00375   dasrt_f = 0;
00376   dasrt_j = 0;
00377   dasrt_cf = 0;
00378 
00379   // Check all the arguments.  Are they the right animals?
00380 
00381   // Here's where I take care of f and j in one shot:
00382 
00383   octave_value f_arg = args(0);
00384 
00385   if (f_arg.is_cell ())
00386     {
00387       Cell c = f_arg.cell_value ();
00388       if (c.length() == 1)
00389         f_arg = c(0);
00390       else if (c.length() == 2)
00391         {
00392           if (c(0).is_function_handle () || c(0).is_inline_function ())
00393             dasrt_f = c(0).function_value ();
00394           else
00395             {
00396               fcn_name = unique_symbol_name ("__dasrt_fcn__");
00397               fname = "function y = ";
00398               fname.append (fcn_name);
00399               fname.append (" (x, xdot, t) y = ");
00400               dasrt_f = extract_function
00401                 (c(0), "dasrt", fcn_name, fname, "; endfunction");
00402             }
00403 
00404           if (dasrt_f)
00405             {
00406               if (c(1).is_function_handle () || c(1).is_inline_function ())
00407                 dasrt_j = c(1).function_value ();
00408               else
00409                 {
00410                   jac_name = unique_symbol_name ("__dasrt_jac__");
00411                   jname = "function jac = ";
00412                   jname.append(jac_name);
00413                   jname.append (" (x, xdot, t, cj) jac = ");
00414                   dasrt_j = extract_function
00415                     (c(1), "dasrt", jac_name, jname, "; endfunction");
00416 
00417                   if (!dasrt_j)
00418                     {
00419                       if (fcn_name.length())
00420                         clear_function (fcn_name);
00421                       dasrt_f = 0;
00422                     }
00423                 }
00424             }
00425         }
00426       else
00427         DASRT_ABORT1 ("incorrect number of elements in cell array");
00428     }
00429 
00430   if (!dasrt_f && ! f_arg.is_cell())
00431     {
00432       if (f_arg.is_function_handle () || f_arg.is_inline_function ())
00433         dasrt_f = f_arg.function_value ();
00434       else
00435         {
00436           switch (f_arg.rows ())
00437             {
00438             case 1:
00439               fcn_name = unique_symbol_name ("__dasrt_fcn__");
00440               fname = "function y = ";
00441               fname.append (fcn_name);
00442               fname.append (" (x, xdot, t) y = ");
00443               dasrt_f = extract_function
00444                 (f_arg, "dasrt", fcn_name, fname, "; endfunction");
00445               break;
00446 
00447             case 2:
00448               {
00449                 string_vector tmp = args(0).all_strings ();
00450 
00451                 if (! error_state)
00452                   {
00453                     fcn_name = unique_symbol_name ("__dasrt_fcn__");
00454                     fname = "function y = ";
00455                     fname.append (fcn_name);
00456                     fname.append (" (x, xdot, t) y = ");
00457                     dasrt_f = extract_function
00458                       (tmp(0), "dasrt", fcn_name, fname, "; endfunction");
00459 
00460                     if (dasrt_f)
00461                       {
00462                         jac_name = unique_symbol_name ("__dasrt_jac__");
00463                         jname = "function jac = ";
00464                         jname.append(jac_name);
00465                         jname.append (" (x, xdot, t, cj) jac = ");
00466                         dasrt_j = extract_function
00467                           (tmp(1), "dasrt", jac_name, jname, "; endfunction");
00468 
00469                         if (! dasrt_j)
00470                           dasrt_f = 0;
00471                       }
00472                   }
00473               }
00474               break;
00475 
00476             default:
00477               DASRT_ABORT1
00478                 ("first arg should be a string or 2-element string array");
00479             }
00480         }
00481     }
00482 
00483   if (error_state || (! dasrt_f))
00484     DASRT_ABORT;
00485 
00486   DAERTFunc func (dasrt_user_f);
00487 
00488   argp++;
00489 
00490   if (args(1).is_function_handle() || args(1).is_inline_function())
00491     {
00492       dasrt_cf = args(1).function_value();
00493 
00494       if (! dasrt_cf)
00495         DASRT_ABORT1 ("expecting function name as argument 2");
00496 
00497       argp++;
00498 
00499       func.set_constraint_function (dasrt_user_cf);
00500     }
00501   else if (args(1).is_string ())
00502     {
00503       dasrt_cf = is_valid_function (args(1), "dasrt", true);
00504       if (! dasrt_cf)
00505         DASRT_ABORT1 ("expecting function name as argument 2");
00506 
00507       argp++;
00508 
00509       func.set_constraint_function (dasrt_user_cf);
00510     }
00511 
00512   ColumnVector state (args(argp++).vector_value ());
00513 
00514   if (error_state)
00515     DASRT_ABORT2 ("expecting state vector as argument %d", argp);
00516 
00517   ColumnVector stateprime (args(argp++).vector_value ());
00518 
00519   if (error_state)
00520     DASRT_ABORT2
00521        ("expecting time derivative of state vector as argument %d", argp);
00522 
00523   ColumnVector out_times (args(argp++).vector_value ());
00524 
00525   if (error_state)
00526     DASRT_ABORT2
00527         ("expecting output time vector as %s argument %d", argp);
00528 
00529   double tzero = out_times (0);
00530 
00531   ColumnVector crit_times;
00532 
00533   bool crit_times_set = false;
00534 
00535   if (argp < nargin)
00536     {
00537       crit_times = ColumnVector (args(argp++).vector_value ());
00538 
00539       if (error_state)
00540         DASRT_ABORT2
00541           ("expecting critical time vector as argument %d", argp);
00542 
00543       crit_times_set = true;
00544     }
00545 
00546   if (dasrt_j)
00547     func.set_jacobian_function (dasrt_user_j);
00548 
00549   DASRT_result output;
00550 
00551   DASRT dae = DASRT (state, stateprime, tzero, func);
00552 
00553   dae.set_options (dasrt_opts);
00554 
00555   if (crit_times_set)
00556     output = dae.integrate (out_times, crit_times);
00557   else
00558     output = dae.integrate (out_times);
00559 
00560   if (fcn_name.length())
00561     clear_function (fcn_name);
00562   if (jac_name.length())
00563     clear_function (jac_name);
00564 
00565   if (! error_state)
00566     {
00567       std::string msg = dae.error_message ();
00568 
00569       retval(4) = msg;
00570       retval(3) = static_cast<double> (dae.integration_state ());
00571 
00572       if (dae.integration_ok ())
00573         {
00574           retval(2) = output.times ();
00575           retval(1) = output.deriv ();
00576           retval(0) = output.state ();
00577         }
00578       else
00579         {
00580           retval(2) = Matrix ();
00581           retval(1) = Matrix ();
00582           retval(0) = Matrix ();
00583 
00584           if (nargout < 4)
00585             error ("dasrt: %s", msg.c_str ());
00586         }
00587     }
00588 
00589   return retval;
00590 }
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines