solution.h

矩阵计算库

//$$ solution.h                      // solve routines

#include "boolean.h"
#include "myexcept.h"

// Solve the equation f(x)=y for x where f is a monotone continuous
// function of x
// Essentially Brent's method

// You need to derive a class from R1_R1 and override "operator()"
// with the function you want to solve.
// Use an object from this class in OneDimSolve

class R1_R1
{
   // the prototype for a Real function of a Real variable
   // you need to derive your function from this one and put in your
   // function for operator() at least. You probably also want to set up a
   // constructor to put in additional parameter values (that won't
   // vary during the solve)

protected:
   Real x;                             // Current x value

public:
   Boolean IsValid() { return TRUE; }  // is the current x value OK
                                       // not used at present
   virtual Real operator()() = 0;      // function value at current x
   virtual void Set(Real X) { x = X; } // set current x
   Boolean IsValid(Real X) { Set(X); return IsValid(); }
                                       // set x, check OK
                                       // not used at present
   Real operator()(Real X) { Set(X); return operator()(); }
                                       // set x, return value
};

class SolutionException : public Exception
{
   static int action;
protected:
   SolutionException() {}
public:
   static long st_type() { return 5; }
   long type() const { return 5; }
   SolutionException(char*);
   SolutionException(const SolutionException&) {}
   static void SetAction(int a) { action=a; }
};

class OneDimSolve
{
   R1_R1& function;                     // reference to the function
   Real acc;                            // accuracy
   int lim;                             // maximum number of iterations

public:
   OneDimSolve(R1_R1& f, Real Acc = 0.0001) : function(f), acc(Acc) {}
                       // f is an R1_R1 function
   Real Solve(Real Y, Real X, Real Dev, int Lim=100);
                       // Solve for x in Y=f(x)
                       // X is the initial trial value of x
                       // X+Dev is the second trial value
                       // program returns a value of x such that
                       // |Y-f(x)| < acc

private:
   Real x[3], y[3];                         // Trial values of X and Y
   int L,C,U,Last;                          // Locations of trial values
   int vpol, hpol;                          // polarities
   Real YY;                                 // target value
   int i;
   void LookAt(int);                        // get new value of function
   Boolean Finish;                          // TRUE if LookAt finds conv.
   void VFlip();
   void HFlip();
   void Flip();
   void State(int I, int J, int K) { L=I; C=J; U=K; }
   void Linear(int, int, int);
   void Quadratic(int, int, int);
};