Mercurial > hg > octave-max
changeset 2850:b7f43611d1e8
[project @ 1997-03-28 21:36:35 by jwe]
author | jwe |
---|---|
date | Fri, 28 Mar 1997 21:39:35 +0000 |
parents | 5338beb20eb9 |
children | b960bd6cbfdf |
files | liboctave/ChangeLog liboctave/LSODE.cc liboctave/LSODE.h src/ChangeLog src/lsode.cc |
diffstat | 5 files changed, 192 insertions(+), 48 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,13 @@ +Fri Mar 28 15:37:09 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * LSODE.h (x_step_limit): New field. + (LSODE_options::init): Initialize it. + (LSODE_options::copy): Copy it. + (LSODE_options::set_step_limit, LSODE_options::step_limit): + New functions. + (LSODE::working_too_hard): Delete. + * LSODE.cc (LSODE::do_integrate): Handle step limit. + Wed Mar 26 15:31:57 1997 John W. Eaton <jwe@bevo.che.wisc.edu> * MArray-b.cc: Delete.
--- a/liboctave/LSODE.cc +++ b/liboctave/LSODE.cc @@ -239,11 +239,12 @@ rwork.elem (5) = (maximum_step_size () >= 0.0) ? maximum_step_size () : 0.0; rwork.elem (6) = (minimum_step_size () >= 0.0) ? minimum_step_size () : 0.0; + if (step_limit () > 0) + iwork.elem (5) = step_limit (); + int *piwork = iwork.fortran_vec (); double *prwork = rwork.fortran_vec (); - working_too_hard = 0; - again: F77_XFCN (lsode, LSODE, (lsode_f, n, xp, t, tout, itol, rel_tol, @@ -268,12 +269,11 @@ break; case -1: // excess work done on this call (perhaps wrong mf). - working_too_hard++; - if (working_too_hard > 20) + if (step_limit () > 0) { (*current_liboctave_error_handler) ("giving up after more than %d steps attempted in lsode", - iwork.elem (5) * 20); + step_limit ()); integration_error = 1; } else
--- a/liboctave/LSODE.h +++ b/liboctave/LSODE.h @@ -64,6 +64,7 @@ x_maximum_step_size = -1.0; x_minimum_step_size = 0.0; x_relative_tolerance = sqrt_eps; + x_step_limit = 500; } void copy (const LSODE_options& opt) @@ -73,6 +74,7 @@ x_maximum_step_size = opt.x_maximum_step_size; x_minimum_step_size = opt.x_minimum_step_size; x_relative_tolerance = opt.x_relative_tolerance; + x_step_limit = opt.x_step_limit; } void set_default_options (void) { init (); } @@ -92,6 +94,9 @@ void set_relative_tolerance (double val) { x_relative_tolerance = (val > 0.0) ? val : ::sqrt (DBL_EPSILON); } + void set_step_limit (int val) + { x_step_limit = val; } + double absolute_tolerance (void) { return x_absolute_tolerance; } @@ -107,6 +112,9 @@ double relative_tolerance (void) { return x_relative_tolerance; } + int step_limit (void) + { return x_step_limit; } + private: double x_absolute_tolerance; @@ -114,6 +122,8 @@ double x_maximum_step_size; double x_minimum_step_size; double x_relative_tolerance; + + int x_step_limit; }; class @@ -164,7 +174,6 @@ int iopt; int liw; int lrw; - int working_too_hard; int sanity_checked; friend int lsode_f (int *neq, double *t, double *y, double *ydot);
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,14 @@ +Fri Mar 28 15:33:11 1997 John W. Eaton <jwe@bevo.che.wisc.edu> + + * lsode.cc (struct LSODE_OPTIONS): Handle integer options. + (print_lsode_option_list, set_lsode_option, show_lsode_option): Ditto. + (lsode_option_table): Add element for step limit. + (lsode_user_jacobian): New function. + (Flsode): Allow function name arg to be a 2-element string array + specifying the function and jacobian function. + + * variables.cc (get_global_value, set_global_value): New functions. + Wed Mar 26 17:08:27 1997 John W. Eaton <jwe@bevo.che.wisc.edu> Implement static variable declaration:
--- a/src/lsode.cc +++ b/src/lsode.cc @@ -43,6 +43,9 @@ // Global pointer for user defined function required by lsode. static tree_fvc *lsode_fcn; +// Global pointer for optional user defined jacobian function used by lsode. +static tree_fvc *lsode_jac; + static LSODE_options lsode_opts; ColumnVector @@ -55,20 +58,11 @@ octave_value_list args; args(1) = t; - if (nstates > 1) - { - Matrix m (nstates, 1); - for (int i = 0; i < nstates; i++) - m (i, 0) = x (i); - octave_value state (m); - args(0) = state; - } - else - { - double d = x (0); - octave_value state (d); - args(0) = state; - } + Matrix m (nstates, 1); + for (int i = 0; i < nstates; i++) + m (i, 0) = x (i); + octave_value state (m); + args(0) = state; if (lsode_fcn) { @@ -94,6 +88,46 @@ return retval; } +Matrix +lsode_user_jacobian (const ColumnVector& x, double t) +{ + Matrix retval; + + int nstates = x.capacity (); + + octave_value_list args; + args(1) = t; + + Matrix m (nstates, 1); + for (int i = 0; i < nstates; i++) + m (i, 0) = x (i); + octave_value state (m); + args(0) = state; + + if (lsode_jac) + { + octave_value_list tmp = lsode_jac->eval (0, 1, args); + + if (error_state) + { + gripe_user_supplied_eval ("lsode"); + return retval; + } + + if (tmp.length () > 0 && tmp(0).is_defined ()) + { + retval = tmp(0).matrix_value (); + + if (error_state || retval.length () == 0) + gripe_user_supplied_eval ("lsode"); + } + else + gripe_user_supplied_eval ("lsode"); + } + + return retval; +} + DEFUN_DLD (lsode, args, nargout, "lsode (F, X0, T_OUT, T_CRIT)\n\ \n\ @@ -114,12 +148,48 @@ return retval; } - lsode_fcn = extract_function - (args(0), "lsode", "__lsode_fcn__", - "function xdot = __lsode_fcn__ (x, t) xdot = ", - "; endfunction"); + octave_value f_arg = args(0); + + switch (f_arg.rows ()) + { + case 1: + lsode_fcn = extract_function + (args(0), "lsode", "__lsode_fcn__", + "function xdot = __lsode_fcn__ (x, t) xdot = ", + "; endfunction"); + break; + + case 2: + { + string_vector tmp = args(0).all_strings (); - if (! lsode_fcn) + if (! error_state) + { + lsode_fcn = extract_function + (tmp(0), "lsode", "__lsode_fcn__", + "function xdot = __lsode_fcn__ (x, t) xdot = ", + "; endfunction"); + + if (lsode_fcn) + { + lsode_jac = extract_function + (tmp(1), "lsode", "__lsode_jac__", + "function jac = __lsode_jac__ (x, t) jac = ", + "; endfunction"); + + if (! lsode_jac) + lsode_fcn = 0; + } + } + } + break; + + default: + error ("lsode: second arg should be a string or 2-element string array"); + break; + } + + if (error_state || ! lsode_fcn) return retval; ColumnVector state = args(1).vector_value (); @@ -158,7 +228,11 @@ int nsteps = out_times.capacity (); ODEFunc func (lsode_user_function); + if (lsode_jac) + func.set_jacobian_function (lsode_user_jacobian); + LSODE ode (state, tzero, func); + ode.copy (lsode_opts); int nstates = state.capacity (); @@ -179,7 +253,9 @@ } typedef void (LSODE_options::*d_set_opt_mf) (double); +typedef void (LSODE_options::*i_set_opt_mf) (int); typedef double (LSODE_options::*d_get_opt_mf) (void); +typedef int (LSODE_options::*i_get_opt_mf) (void); #define MAX_TOKENS 3 @@ -190,7 +266,9 @@ int min_len[MAX_TOKENS + 1]; int min_toks_to_match; d_set_opt_mf d_set_fcn; + i_set_opt_mf i_set_fcn; d_get_opt_mf d_get_fcn; + i_get_opt_mf i_get_fcn; }; static LSODE_OPTIONS lsode_option_table [] = @@ -198,37 +276,43 @@ { "absolute tolerance", { "absolute", "tolerance", 0, 0, }, { 1, 0, 0, 0, }, 1, - LSODE_options::set_absolute_tolerance, - LSODE_options::absolute_tolerance, }, + LSODE_options::set_absolute_tolerance, 0, + LSODE_options::absolute_tolerance, 0, }, { "initial step size", { "initial", "step", "size", 0, }, { 1, 0, 0, 0, }, 1, - LSODE_options::set_initial_step_size, - LSODE_options::initial_step_size, }, + LSODE_options::set_initial_step_size, 0, + LSODE_options::initial_step_size, 0, }, { "maximum step size", { "maximum", "step", "size", 0, }, { 2, 0, 0, 0, }, 1, - LSODE_options::set_maximum_step_size, - LSODE_options::maximum_step_size, }, + LSODE_options::set_maximum_step_size, 0, + LSODE_options::maximum_step_size, 0, }, { "minimum step size", { "minimum", "step", "size", 0, }, { 2, 0, 0, 0, }, 1, - LSODE_options::set_minimum_step_size, - LSODE_options::minimum_step_size, }, + LSODE_options::set_minimum_step_size, 0, + LSODE_options::minimum_step_size, 0, }, { "relative tolerance", { "relative", "tolerance", 0, 0, }, { 1, 0, 0, 0, }, 1, - LSODE_options::set_relative_tolerance, - LSODE_options::relative_tolerance, }, + LSODE_options::set_relative_tolerance, 0, + LSODE_options::relative_tolerance, 0, }, + + { "step limit", + { "step", "limit", 0, 0, }, + { 1, 0, 0, 0, }, 1, + 0, LSODE_options::set_step_limit, + 0, LSODE_options::step_limit, }, { 0, { 0, 0, 0, 0, }, { 0, 0, 0, 0, }, 0, - 0, 0, }, + 0, 0, 0, 0, }, }; static void @@ -247,13 +331,22 @@ while ((keyword = list->keyword) != 0) { os.form (" %-40s ", keyword); - - double val = (lsode_opts.*list->d_get_fcn) (); - if (val < 0.0) - os << "computed automatically"; + if (list->d_get_fcn) + { + double val = (lsode_opts.*list->d_get_fcn) (); + if (val < 0.0) + os << "computed automatically"; + else + os << val; + } else - os << val; - + { + int val = (lsode_opts.*list->i_get_fcn) (); + if (val < 0) + os << "infinite"; + else + os << val; + } os << "\n"; list++; } @@ -271,8 +364,18 @@ if (keyword_almost_match (list->kw_tok, list->min_len, keyword, list->min_toks_to_match, MAX_TOKENS)) { - (lsode_opts.*list->d_set_fcn) (val); - + if (list->d_set_fcn) + (lsode_opts.*list->d_set_fcn) (val); + else + { + if (xisnan (val)) + { + error ("lsode_options: %s: expecting integer, found NaN", + keyword.c_str ()); + } + else + (lsode_opts.*list->i_set_fcn) (NINT (val)); + } return; } list++; @@ -293,11 +396,22 @@ if (keyword_almost_match (list->kw_tok, list->min_len, keyword, list->min_toks_to_match, MAX_TOKENS)) { - double val = (lsode_opts.*list->d_get_fcn) (); - if (val < 0.0) - retval = "computed automatically"; + if (list->d_get_fcn) + { + double val = (lsode_opts.*list->d_get_fcn) (); + if (val < 0.0) + retval = "computed automatically"; + else + retval = val; + } else - retval = val; + { + int val = (lsode_opts.*list->i_get_fcn) (); + if (val < 0) + retval = "infinite"; + else + retval = static_cast<double> (val); + } return retval; }