adams_pece_decl.h

时间序列分析的信号处理MATLAB工具包

#ifndef ADAM_PECE_IMPL_H
#define ADAM_PECE_IMPL_H

#include <math.h>

class NumericalAlgorithm {
	public:
		static double machine_precision(); 	// return machine precision		
};

#ifdef FORTRAN_LIBRARY_AVAILABLE
// If the native fortran library is available, use the mathmatical
// functions for total compatibility 
// In this case, you have to link with -lftn or -lF77 and -lI77
extern "C" {		
	long i_sign(long*, long *);
	double d_sign(double*, double*);
	double pow_dd(double*, double*);
}
#else
// Otherwise, give substitutes for the original functions
class f2clib {	
	public:
		inline static double d_sign(double *a, double *b)
		{
			const double x = (*a >= 0 ? *a : - *a);
			return( *b >= 0 ? x : -x);
		}
		inline static long i_sign(long *a, long *b)
		{               
			const long x = (*a >= 0 ? *a : - *a);
			return( *b >= 0 ? x : -x);
		}
		inline static double pow_dd(double *ap, double *bp)
		{
			return(pow(*ap, *bp) );
		}
		template<class T> static T abs(const T x) { return (x >= 0) ? x : -x; }
		static double dabs(const double x) { return (double) abs(x); }	
		template<class T> static T min(const T a, const T b) { return ((a) <= (b) ? (a) : (b)); }
		template<class T> static T max(const T a, const T b) { return ((a) >= (b) ? (a) : (b)); }
		static double dmin(const double a, const double b) { return (double) min(a,b); }
		static double dmax(const double a, const double b) { return (double) max(a,b); }		
};
#endif

template<class ExplicitODE>
#ifdef FORTRAN_LIBRARY_AVAILABLE
class Adams : private NumericalAlgorithm {
#else
class Adams : public NumericalAlgorithm, private f2clib {
#endif
	private:
		const ExplicitODE& f;
		const long neqn;
		
		// machine precision parameters
		const double u; 
		const double twou;
		const double fouru;
		
		// constants generated by f2c
		long c__4;
		long c__1;
		double c_b11;
		
		long iwork_c[5];
		double* work_c; 
		
		int de(const long neqn, double *y, double *t, double *tout, double *relerr, double *abserr, long *iflag, double *yy, double *wt, 
			double *p, double *yp, double *ypout, double *phi, double *alpha, double *beta, double *sig, double *v, double *w, 
			double *g, long *phase1, double *psi, double *x, double *h__, double *hold, long *start, double *told, double *delsgn, 
			long *ns, long *nornd, long *k, long *kold, long *isnold);
			
		int step(double *x, double *y, const long neqn, double *h__, double *eps, double *wt, long *start, double *hold, long *k, 
			long *kold, long *crash, double *phi, double *p, double *yp, double *psi, double *alpha, double *beta, double *sig, 
			double *v, double *w, double *g, long *phase1, long *ns, long *nornd);
		
		int intrp(const double x, const double *y, const double xout, double *yout, double *ypout, const long neqn, const long kold, double *phi, double *psi);
			
	public:
		Adams(const ExplicitODE& system) : f(system), neqn(system.get_number_equations()), u(machine_precision()), twou(2*u), fouru(4*u), 
		c__4(4), c__1(1),  c_b11(1.0), work_c(0) {
			work_c = new double [100+21*neqn];
		};
		~Adams() {
			delete[] work_c;
		};
		
		int ode(double *t, double *tout, double *relerr, double *abserr, long *iflag);	
		
//    double precision subroutine ode integrates a system of neqn
//    first order ordinary differential equations of the form:
//              dy(i)/dt = f(t,y(1),y(2),...,y(neqn))
//              y(i) given at  t .
//    the subroutine integrates from  t  to  tout .  on return the
//    parameters in the call list are set for continuing the integration.
//    the user has only to define a new value  tout  and call  ode  again.
//
//    the differential equations are actually solved by a suite of codes
//    de ,  step , and  intrp .  ode calls on the routines  step  and  intrp
//    to advance the integration and to interpolate at output points.
//    step  uses a modified divided difference form of the adams pece
//    formulas and local extrapolation.  it adjusts the order and step
//    size to control the local error per unit step in a generalized
//    sense.  normally each call to  step  advances the solution one step
//    in the direction of  tout .  for reasons of efficiency  de
//    integrates beyond  tout  internally, though never beyond
//    t+10*(tout-t), and calls  intrp  to interpolate the solution at
//    tout .  an option is provided to stop the integration at  tout  but
//    it should be used only if it is impossible to continue the
//    integration beyond  tout .
//
//    this code is completely explained and documented in the text,
//    computer solution of ordinary differential equations:  the initial
//    value problem  by l. f. shampine and m. k. gordon.
//
//    the parameters represent:
//       f -- double precision subroutine f(t,y,yp) to evaluate
//                 derivatives yp(i)=dy(i)/dt
//       neqn -- number of equations to be integrated (integer*4)
//       y(*) -- solution vector at t                 (real*8)
//       t -- independent variable                    (real*8)
//       tout -- point at which solution is desired   (real*8)
//       relerr,abserr -- relative and absolute error tolerances for local
//            error test (real*8).  at each step the code requires
//              dabs(local error) < dabs(y)*relerr + abserr
//            for each component of the local error and solution vectors
//       iflag -- indicates status of integration     (integer*4)
//
//    first call to ode --
//
//
//       t -- starting point of integration
//       tout -- point at which solution is desired
//       relerr,abserr -- relative and absolute local error tolerances
//       iflag -- +1,-1.  indicator to initialize the code.  normal input
//            is +1.  the user should set iflag=-1 only if it is
//            impossible to continue the integration beyond  tout .
//    all parameters except tout  may be altered by the
//    code on output so must be variables in the calling program.
//
//    output from  ode  --
//
//       t -- last point reached in integration.  normal return has
//            t = tout .
//       tout -- unchanged
//       relerr,abserr -- normal return has tolerances unchanged.  iflag=3
//            signals tolerances increased
//       iflag = 2 -- normal return.  integration reached  tout
//             = 3 -- integration did not reach  tout  because error
//                    tolerances too small.  relerr ,  abserr  increased
//                    appropriately for continuing
//             = 4 -- integration did not reach  tout  because more than
//                    1000 steps needed
//             = 5 -- integration did not reach  tout  because equations
//                    appear to be stiff
//             = 6 -- invalid input parameters (fatal error)
//            the value of  iflag  is returned negative when the input
//            value is negative and the integration does not reach  tout ,
//            i.e., -3, -4, -5.
//
//    subsequent calls to  ode --
//
//    subroutine  ode  returns with all information needed to continue
//    the integration.  if the integration reached  tout , the user need
//    only define a new  tout  and call again.  if the integration did not
//    reach  tout  and the user wants to continue, he just calls again.
//    the output value of  iflag  is the appropriate input value for
//    subsequent calls.  the only situation in which it should be altered
//    is to stop the integration internally at the new  tout , i.e.,
//    change output  iflag=2  to input  iflag=-2 .  error tolerances may
//    be changed by the user before continuing.  all other parameters must
//    remain unchanged.

			
};

	

#endif