81 if (tmp.
length () > 0 && tmp(0).is_defined ())
85 warning (
"lsode: ignoring imaginary part returned from user-supplied function");
120 if (tmp.
length () > 0 && tmp(0).is_defined ())
124 warning (
"lsode: ignoring imaginary part returned from user-supplied jacobian function");
128 retval = tmp(0).matrix_value ();
140 #define LSODE_ABORT() \
143 #define LSODE_ABORT1(msg) \
146 ::error ("lsode: " msg); \
151 #define LSODE_ABORT2(fmt, arg) \
154 ::error ("lsode: " fmt, arg); \
159 DEFUN (lsode, args, nargout,
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\
165 The set of differential equations to solve is\n\
167 $$ {dx \\over dt} = f (x, t) $$\n\
169 $$ x(t_0) = x_0 $$\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\
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\
200 @var{xdot} = f (@var{x}, @var{t})\n\
204 in which @var{xdot} and @var{x} are vectors and @var{t} is a scalar.\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\
213 @var{jac} = j (@var{x}, @var{t})\n\
217 in which @var{jac} is the matrix of partial derivatives\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\
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\
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\
232 & {\\partial f_3 \\over \\partial x_N} \\cr}\\right]$$\n\
238 | df_1 df_1 df_1 |\n\
239 | ---- ---- ... ---- |\n\
240 | dx_1 dx_2 dx_N |\n\
242 | df_2 df_2 df_2 |\n\
243 | ---- ---- ... ---- |\n\
244 df_i | dx_1 dx_2 dx_N |\n\
250 | df_N df_N df_N |\n\
251 | ---- ---- ... ---- |\n\
252 | dx_1 dx_2 dx_N |\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\
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\
266 After a successful computation, the value of @var{istate} will be 2\n\
267 (consistent with the Fortran version of @sc{lsode}).\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\
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\
290 int nargin = args.length ();
292 if (nargin > 2 && nargin < 5 && nargout < 4)
294 std::string fcn_name, fname, jac_name, jname;
305 else if (c.
length () == 2)
307 if (c(0).is_function_handle () || c(0).is_inline_function ())
308 lsode_fcn = c(0).function_value ();
312 fname =
"function y = ";
313 fname.append (fcn_name);
314 fname.append (
" (x, t) y = ");
321 if (c(1).is_function_handle () || c(1).is_inline_function ())
322 lsode_jac = c(1).function_value ();
326 jname =
"function jac = ";
327 jname.append (jac_name);
328 jname.append (
" (x, t) jac = ");
330 jname,
"; endfunction");
334 if (fcn_name.length ())
342 LSODE_ABORT1 (
"incorrect number of elements in cell array");
345 if (!lsode_fcn && ! f_arg.
is_cell ())
351 switch (f_arg.
rows ())
357 fname =
"function y = ";
358 fname.append (fcn_name);
359 fname.append (
" (x, t) y = ");
361 fname,
"; endfunction");
373 fname =
"function y = ";
374 fname.append (fcn_name);
375 fname.append (
" (x, t) y = ");
377 fname,
"; endfunction");
382 jname =
"function jac = ";
383 jname.append (jac_name);
384 jname.append (
" (x, t) jac = ");
391 if (fcn_name.length ())
402 (
"first arg should be a string or 2-element string array");
413 LSODE_ABORT1 (
"expecting state vector as second argument");
418 LSODE_ABORT1 (
"expecting output time vector as third argument");
422 int crit_times_set = 0;
428 LSODE_ABORT1 (
"expecting critical time vector as fourth argument");
433 double tzero = out_times (0);
439 LSODE ode (state, tzero, func);
445 output = ode.
integrate (out_times, crit_times);
449 if (fcn_name.length ())
451 if (jac_name.length ())
468 error (
"lsode: %s", msg.c_str ());
static bool warned_jac_imaginary
virtual octave_value_list do_multi_index_op(int nargout, const octave_value_list &idx)
octave_idx_type rows(void) const
bool integration_ok(void) const
static uint32_t state[624]
OCTINTERP_API void print_usage(void)
octave_idx_type length(void) const
Matrix lsode_user_jacobian(const ColumnVector &x, double t)
virtual ColumnVector integrate(double tt)
static bool warned_fcn_imaginary
#define DEFUN(name, args_name, nargout_name, doc)
void error(const char *fmt,...)
ColumnVector lsode_user_function(const ColumnVector &x, double t)
void set_options(const LSODE_options &opt)
bool is_function_handle(void) const
string_vector all_strings(bool pad=false) const
std::string unique_symbol_name(const std::string &basename)
Cell cell_value(void) const
void clear_function(const std::string &nm)
ODEFunc & set_jacobian_function(ODEJacFunc j)
static octave_function * lsode_jac
bool is_inline_function(void) const
#define LSODE_ABORT1(msg)
octave_function * function_value(bool silent=false) const
void warning(const char *fmt,...)
static octave_function * lsode_fcn
octave_idx_type length(void) const
Number of elements in the array.
static LSODE_options lsode_opts
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)
octave_idx_type integration_state(void) const
void gripe_user_supplied_eval(const char *name)
F77_RET_T const double * x
std::string error_message(void) const