Mercurial > hg > octave-lyh
changeset 1842:0574a1f3a273
[project @ 1996-02-03 11:44:02 by jwe]
author | jwe |
---|---|
date | Sat, 03 Feb 1996 11:44:02 +0000 |
parents | fc5667a20dd2 |
children | 88c5728ae7ae |
files | liboctave/DAE.h liboctave/DASSL.cc liboctave/LSODE.cc liboctave/ODE.h |
diffstat | 4 files changed, 251 insertions(+), 345 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/DAE.h +++ b/liboctave/DAE.h @@ -1,7 +1,7 @@ // DAE.h -*- C++ -*- /* -Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton +Copyright (C) 1996 John W. Eaton This file is part of Octave. @@ -28,61 +28,49 @@ #pragma interface #endif -#include "dColVector.h" -#include "ODE.h" #include "DAEFunc.h" +#include "base-de.h" -class DAE : public ODE, public DAEFunc +class DAE : public base_diff_eqn, public DAEFunc { public: - DAE (void); + DAE (void) + : base_diff_eqn (), DAEFunc (), xdot () { } - DAE (int); + DAE (const ColumnVector& x, double t, DAEFunc& f) + : base_diff_eqn (x, t), DAEFunc (f), xdot (x.capacity (), 0.0) { } - DAE (const ColumnVector& x, double time, DAEFunc& f); + DAE (const ColumnVector& x, const ColumnVector& xxdot, + double t, DAEFunc& f); - DAE (const ColumnVector& x, const ColumnVector& xdot, - double time, DAEFunc& f); - - ~DAE (void); + DAE (const DAE& a) + : base_diff_eqn (a), DAEFunc (a), xdot (a.xdot) { } - ColumnVector deriv (void); - - virtual void initialize (const ColumnVector& x, double t); - virtual void initialize (const ColumnVector& x, - const ColumnVector& xdot, double t); + DAE& operator = (const DAE& a) + { + if (this != &a) + { + base_diff_eqn::operator = (a); + DAEFunc::operator = (a); - ColumnVector integrate (double t); + xdot = a.xdot; + } + return *this; + } + + ~DAE (void) { } - Matrix integrate (const ColumnVector& tout, Matrix& xdot_out); - Matrix integrate (const ColumnVector& tout, Matrix& xdot_out, - const ColumnVector& tcrit); + ColumnVector state_derivative (void) { return xdot; } + + void initialize (const ColumnVector& x, double t); + + void initialize (const ColumnVector& x, const ColumnVector& xxdot, + double t); protected: - // Some of this is probably too closely related to DASSL, but hey, - // this is just a first attempt... - ColumnVector xdot; - -private: - - int integration_error; - int restart; - int liw; - int lrw; - int idid; - int *info; - int *iwork; - double *rwork; - - friend int ddassl_j (double *time, double *state, double *deriv, - double *pd, double *cj, double *rpar, int *ipar); - - friend int ddassl_f (double *time, double *state, double *deriv, - double *delta, int *ires, double *rpar, int *ipar); - }; #endif
--- a/liboctave/DASSL.cc +++ b/liboctave/DASSL.cc @@ -1,7 +1,7 @@ -// DAE.cc -*- C++ -*- +// DASSL.cc -*- C++ -*- /* -Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton +Copyright (C) 1996 John W. Eaton This file is part of Octave. @@ -29,7 +29,10 @@ #include <config.h> #endif -#include "DAE.h" +#include <cfloat> +#include <cmath> + +#include "DASSL.h" #include "f77-uscore.h" #include "lo-error.h" @@ -52,20 +55,11 @@ static DAEFunc::DAEJacFunc user_jac; static int nn; -DAE::DAE (void) +DASSL::DASSL (void) : DAE () { - n = 0; - t = 0.0; - stop_time_set = 0; stop_time = 0.0; - integration_error = 0; - restart = 1; - - DAEFunc::set_function (0); - DAEFunc::set_jacobian_function (0); - liw = 0; lrw = 0; @@ -77,20 +71,14 @@ info [i] = 0; } -DAE::DAE (int size) +DASSL::DASSL (const ColumnVector& state, double time, DAEFunc& f) + : DAE (state, time, f) { - n = size; - t = 0.0; + n = size (); stop_time_set = 0; stop_time = 0.0; - integration_error = 0; - restart = 1; - - DAEFunc::set_function (0); - DAEFunc::set_jacobian_function (0); - liw = 20 + n; lrw = 40 + 9*n + n*n; @@ -102,49 +90,11 @@ info [i] = 0; } -DAE::DAE (const ColumnVector& state, double time, DAEFunc& f) +DASSL::DASSL (const ColumnVector& state, const ColumnVector& deriv, + double time, DAEFunc& f) + : DAE (state, deriv, time, f) { - n = state.capacity (); - t = time; - x = state; - xdot.resize (n, 0.0); - - stop_time_set = 0; - stop_time = 0.0; - - integration_error = 0; - restart = 1; - - DAEFunc::set_function (f.function ()); - DAEFunc::set_jacobian_function (f.jacobian_function ()); - - liw = 20 + n; - lrw = 40 + 9*n + n*n; - - info = new int [15]; - iwork = new int [liw]; - rwork = new double [lrw]; - - for (int i = 0; i < 15; i++) - info [i] = 0; -} - -DAE::DAE (const ColumnVector& state, const ColumnVector& deriv, - double time, DAEFunc& f) -{ - if (deriv.capacity () != state.capacity ()) - { - (*current_liboctave_error_handler) - ("x, xdot size mismatch in DAE constructor"); - n = 0; - t = 0.0; - return; - } - - n = state.capacity (); - t = time; - xdot = deriv; - x = state; + n = size (); stop_time_set = 0; stop_time = 0.0; @@ -163,39 +113,31 @@ info [i] = 0; } -DAE::~DAE (void) +DASSL::~DASSL (void) { delete info; delete rwork; delete iwork; } -ColumnVector -DAE::deriv (void) +void +DASSL::force_restart (void) { - return xdot; + restart = 1; + integration_error = 0; } void -DAE::initialize (const ColumnVector& state, double time) +DASSL::set_stop_time (double t) { - integration_error = 0; - restart = 1; - x = state; - int nx = x.capacity (); - xdot.resize (nx, 0.0); - t = time; + stop_time_set = 1; + stop_time = t; } void -DAE::initialize (const ColumnVector& state, - const ColumnVector& deriv, double time) +DASSL::clear_stop_time (void) { - integration_error = 0; - restart = 1; - xdot = deriv; - x = state; - t = time; + stop_time_set = 0; } int @@ -255,7 +197,7 @@ } ColumnVector -DAE::integrate (double tout) +DASSL::do_integrate (double tout) { integration_error = 0; @@ -363,7 +305,14 @@ } Matrix -DAE::integrate (const ColumnVector& tout, Matrix& xdot_out) +DASSL::do_integrate (const ColumnVector& tout) +{ + Matrix dummy; + return integrate (tout, dummy); +} + +Matrix +DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out) { Matrix retval; int n_out = tout.capacity (); @@ -381,7 +330,7 @@ for (int j = 1; j < n_out; j++) { - ColumnVector x_next = integrate (tout.elem (j)); + ColumnVector x_next = do_integrate (tout.elem (j)); if (integration_error) return retval; @@ -398,8 +347,8 @@ } Matrix -DAE::integrate (const ColumnVector& tout, Matrix& xdot_out, - const ColumnVector& tcrit) +DASSL::integrate (const ColumnVector& tout, Matrix& xdot_out, + const ColumnVector& tcrit) { Matrix retval; int n_out = tout.capacity (); @@ -469,7 +418,7 @@ i_out++; } - ColumnVector x_next = integrate (t_out); + ColumnVector x_next = do_integrate (t_out); if (integration_error) return retval; @@ -499,6 +448,116 @@ return retval; } +DASSL_options::DASSL_options (void) +{ + init (); +} + +DASSL_options::DASSL_options (const DASSL_options& opt) +{ + copy (opt); +} + +DASSL_options& +DASSL_options::operator = (const DASSL_options& opt) +{ + if (this != &opt) + copy (opt); + + return *this; +} + +DASSL_options::~DASSL_options (void) +{ +} + +void +DASSL_options::init (void) +{ + double sqrt_eps = sqrt (DBL_EPSILON); + x_absolute_tolerance = sqrt_eps; + x_initial_step_size = -1.0; + x_maximum_step_size = -1.0; + x_minimum_step_size = 0.0; + x_relative_tolerance = sqrt_eps; +} + +void +DASSL_options::copy (const DASSL_options& opt) +{ + x_absolute_tolerance = opt.x_absolute_tolerance; + x_initial_step_size = opt.x_initial_step_size; + 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; +} + +void +DASSL_options::set_default_options (void) +{ + init (); +} + +void +DASSL_options::set_absolute_tolerance (double val) +{ + x_absolute_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); +} + +void +DASSL_options::set_initial_step_size (double val) +{ + x_initial_step_size = (val >= 0.0) ? val : -1.0; +} + +void +DASSL_options::set_maximum_step_size (double val) +{ + x_maximum_step_size = (val >= 0.0) ? val : -1.0; +} + +void +DASSL_options::set_minimum_step_size (double val) +{ + x_minimum_step_size = (val >= 0.0) ? val : 0.0; +} + +void +DASSL_options::set_relative_tolerance (double val) +{ + x_relative_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); +} + +double +DASSL_options::absolute_tolerance (void) +{ + return x_absolute_tolerance; +} + +double +DASSL_options::initial_step_size (void) +{ + return x_initial_step_size; +} + +double +DASSL_options::maximum_step_size (void) +{ + return x_maximum_step_size; +} + +double +DASSL_options::minimum_step_size (void) +{ + return x_minimum_step_size; +} + +double +DASSL_options::relative_tolerance (void) +{ + return x_relative_tolerance; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/LSODE.cc +++ b/liboctave/LSODE.cc @@ -1,4 +1,4 @@ -// ODE.cc -*- C++ -*- +// LSODE.cc -*- C++ -*- /* Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton @@ -34,7 +34,7 @@ #include <iostream.h> -#include "ODE.h" +#include "LSODE.h" #include "f77-uscore.h" #include "lo-error.h" @@ -55,10 +55,9 @@ static ODEFunc::ODEJacFunc user_jac; static ColumnVector *tmp_x; -ODE::ODE (void) +LSODE::LSODE (void) : ODE (), LSODE_options () { - n = 0; - t = 0.0; + n = size (); stop_time_set = 0; stop_time = 0.0; @@ -81,15 +80,12 @@ iwork[i] = 0; rwork[i] = 0.0; } - - fun = 0; - jac = 0; } -ODE::ODE (int size) +LSODE::LSODE (const ColumnVector& state, double time, const ODEFunc& f) + : ODE (state, time, f), LSODE_options () { - n = size; - t = 0.0; + n = size (); stop_time_set = 0; stop_time = 0.0; @@ -112,49 +108,33 @@ iwork[i] = 0; rwork[i] = 0.0; } - - fun = 0; - jac = 0; } -ODE::ODE (const ColumnVector& state, double time, const ODEFunc& f) -{ - n = state.capacity (); - t = time; - x = state; - - stop_time_set = 0; - stop_time = 0.0; - - integration_error = 0; - restart = 1; - - istate = 1; - itol = 1; - itask = 1; - iopt = 1; - - liw = 20 + n; - lrw = 22 + n * (9 + n); - - iwork = new int [liw]; - rwork = new double [lrw]; - for (int i = 4; i < 9; i++) - { - iwork[i] = 0; - rwork[i] = 0.0; - } - - fun = f.function (); - jac = f.jacobian_function (); -} - -ODE::~ODE (void) +LSODE::~LSODE (void) { delete [] rwork; delete [] iwork; } +void +LSODE::force_restart (void) +{ + restart = 1; +} + +void +LSODE::set_stop_time (double time) +{ + stop_time_set = 1; + stop_time = time; +} + +void +LSODE::clear_stop_time (void) +{ + stop_time_set = 0; +} + int lsode_f (const int& neq, const double& time, double *, double *deriv, int& ierr) @@ -198,7 +178,7 @@ } ColumnVector -ODE::integrate (double tout) +LSODE::do_integrate (double tout) { if (jac) method_flag = 21; @@ -214,8 +194,8 @@ // and copy. tmp_x = &x; - user_fun = fun; - user_jac = jac; + user_fun = function (); + user_jac = jacobian_function (); // Try 5000 steps before giving up. @@ -295,8 +275,9 @@ return x; } +#if 0 void -ODE::integrate (int nsteps, double tstep, ostream& s) +LSODE::integrate (int nsteps, double tstep, ostream& s) { int time_to_quit = 0; double tout = t; @@ -320,9 +301,10 @@ return; } } +#endif Matrix -ODE::integrate (const ColumnVector& tout) +LSODE::do_integrate (const ColumnVector& tout) { Matrix retval; int n_out = tout.capacity (); @@ -336,7 +318,7 @@ for (int j = 1; j < n_out; j++) { - ColumnVector x_next = integrate (tout.elem (j)); + ColumnVector x_next = do_integrate (tout.elem (j)); if (integration_error) return retval; @@ -350,7 +332,7 @@ } Matrix -ODE::integrate (const ColumnVector& tout, const ColumnVector& tcrit) +LSODE::integrate (const ColumnVector& tout, const ColumnVector& tcrit) { Matrix retval; int n_out = tout.capacity (); @@ -416,7 +398,7 @@ i_out++; } - ColumnVector x_next = integrate (t_out); + ColumnVector x_next = do_integrate (t_out); if (integration_error) return retval; @@ -433,7 +415,7 @@ } else { - retval = integrate (tout); + retval = do_integrate (tout); if (integration_error) return retval; @@ -443,62 +425,18 @@ return retval; } -int -ODE::size (void) const -{ - return n; -} - -ColumnVector -ODE::state (void) const -{ - return x; -} - -double ODE::time (void) const -{ - return t; -} - -void -ODE::force_restart (void) -{ - restart = 1; -} - -void -ODE::initialize (const ColumnVector& state, double time) -{ - restart = 1; - x = state; - t = time; -} - -void -ODE::set_stop_time (double time) -{ - stop_time_set = 1; - stop_time = time; -} - -void -ODE::clear_stop_time (void) -{ - stop_time_set = 0; -} - -ODE_options::ODE_options (void) +LSODE_options::LSODE_options (void) { init (); } -ODE_options::ODE_options (const ODE_options& opt) +LSODE_options::LSODE_options (const LSODE_options& opt) { copy (opt); } -ODE_options& -ODE_options::operator = (const ODE_options& opt) +LSODE_options& +LSODE_options::operator = (const LSODE_options& opt) { if (this != &opt) copy (opt); @@ -506,12 +444,12 @@ return *this; } -ODE_options::~ODE_options (void) +LSODE_options::~LSODE_options (void) { } void -ODE_options::init (void) +LSODE_options::init (void) { double sqrt_eps = sqrt (DBL_EPSILON); x_absolute_tolerance = sqrt_eps; @@ -522,7 +460,7 @@ } void -ODE_options::copy (const ODE_options& opt) +LSODE_options::copy (const LSODE_options& opt) { x_absolute_tolerance = opt.x_absolute_tolerance; x_initial_step_size = opt.x_initial_step_size; @@ -532,67 +470,67 @@ } void -ODE_options::set_default_options (void) +LSODE_options::set_default_options (void) { init (); } void -ODE_options::set_absolute_tolerance (double val) +LSODE_options::set_absolute_tolerance (double val) { x_absolute_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); } void -ODE_options::set_initial_step_size (double val) +LSODE_options::set_initial_step_size (double val) { x_initial_step_size = (val >= 0.0) ? val : -1.0; } void -ODE_options::set_maximum_step_size (double val) +LSODE_options::set_maximum_step_size (double val) { x_maximum_step_size = (val >= 0.0) ? val : -1.0; } void -ODE_options::set_minimum_step_size (double val) +LSODE_options::set_minimum_step_size (double val) { x_minimum_step_size = (val >= 0.0) ? val : 0.0; } void -ODE_options::set_relative_tolerance (double val) +LSODE_options::set_relative_tolerance (double val) { x_relative_tolerance = (val > 0.0) ? val : sqrt (DBL_EPSILON); } double -ODE_options::absolute_tolerance (void) +LSODE_options::absolute_tolerance (void) { return x_absolute_tolerance; } double -ODE_options::initial_step_size (void) +LSODE_options::initial_step_size (void) { return x_initial_step_size; } double -ODE_options::maximum_step_size (void) +LSODE_options::maximum_step_size (void) { return x_maximum_step_size; } double -ODE_options::minimum_step_size (void) +LSODE_options::minimum_step_size (void) { return x_minimum_step_size; } double -ODE_options::relative_tolerance (void) +LSODE_options::relative_tolerance (void) { return x_relative_tolerance; }
--- a/liboctave/ODE.h +++ b/liboctave/ODE.h @@ -1,7 +1,7 @@ // ODE.h -*- C++ -*- /* -Copyright (C) 1992, 1993, 1994, 1995 John W. Eaton +Copyright (C) 1996 John W. Eaton This file is part of Octave. @@ -24,112 +24,33 @@ #if !defined (octave_ODE_h) #define octave_ODE_h 1 -#if defined (__GNUG__) -#pragma interface -#endif - -class ostream; - -#include "dMatrix.h" -#include "dColVector.h" #include "ODEFunc.h" - -class ODE_options -{ - public: - - ODE_options (void); - ODE_options (const ODE_options& opt); - - ODE_options& operator = (const ODE_options& opt); - - ~ODE_options (void); - - void init (void); - void copy (const ODE_options& opt); +#include "base-de.h" - void set_default_options (void); - - void set_absolute_tolerance (double); - void set_initial_step_size (double); - void set_maximum_step_size (double); - void set_minimum_step_size (double); - void set_relative_tolerance (double); - - double absolute_tolerance (void); - double initial_step_size (void); - double maximum_step_size (void); - double minimum_step_size (void); - double relative_tolerance (void); - - private: - - double x_absolute_tolerance; - double x_initial_step_size; - double x_maximum_step_size; - double x_minimum_step_size; - double x_relative_tolerance; -}; - -class ODE : public ODEFunc, public ODE_options +class ODE : public base_diff_eqn, public ODEFunc { public: - ODE (void); - - ODE (int n); - - ODE (const ColumnVector& state, double time, const ODEFunc& f); - - virtual ~ODE (void); - - virtual int size (void) const; - virtual ColumnVector state (void) const; - virtual double time (void) const; + ODE (void) + : base_diff_eqn (), ODEFunc () { } - virtual void force_restart (void); - virtual void initialize (const ColumnVector& x, double t); - virtual void set_stop_time (double t); - virtual void clear_stop_time (void); - - virtual ColumnVector integrate (double t); + ODE (const ColumnVector& state, double time, const ODEFunc& f) + : base_diff_eqn (state, time), ODEFunc (f) { } - void integrate (int nsteps, double tstep, ostream& s); - - Matrix integrate (const ColumnVector& tout); - Matrix integrate (const ColumnVector& tout, const ColumnVector& tcrit); - -protected: + ODE (const ODE& a) + : base_diff_eqn (a), ODEFunc (a) { } - // Some of this is probably too closely related to LSODE, but hey, - // this is just a first attempt... - - int n; - double t; - ColumnVector x; - - double stop_time; - int stop_time_set; - -private: + ODE& operator = (const ODE& a) + { + if (this != &a) + { + base_diff_eqn::operator = (a); + ODEFunc::operator = (a); + } + return *this; + } - int integration_error; - int restart; - int method_flag; - int *iwork; - double *rwork; - int istate; - int itol; - int itask; - int iopt; - int liw; - int lrw; - - friend int lsode_f (int *neq, double *t, double *y, double *ydot); - - friend int lsode_j (int *neq, double *t, double *y, int *ml, int *mu, - double *pd, int *nrowpd); - + ~ODE (void) { } }; #endif