Mercurial > hg > octave-jordi

changeset 1842:0574a1f3a273

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[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