$search
00001 00002 00003 00006 00007 #include "myexcept.h" 00008 00009 #ifdef use_namespace 00010 namespace RBD_COMMON { 00011 #endif 00012 00013 00014 // Solve the equation f(x)=y for x where f is a monotone continuous 00015 // function of x 00016 // Essentially Brent s method 00017 00018 // You need to derive a class from R1_R1 and override "operator()" 00019 // with the function you want to solve. 00020 // Use an object from this class in OneDimSolve 00021 00022 class R1_R1 00023 { 00024 // the prototype for a Real function of a Real variable 00025 // you need to derive your function from this one and put in your 00026 // function for operator() at least. You probably also want to set up a 00027 // constructor to put in additional parameter values (e.g. that will not 00028 // vary during a solve) 00029 00030 protected: 00031 Real x; // Current x value 00032 bool xSet; // true if a value assigned to x 00033 00034 public: 00035 Real minX, maxX; // range of value x 00036 bool minXinf, maxXinf; // true if these are infinite 00037 R1_R1() : xSet(false), minXinf(true), maxXinf(true) {} 00038 virtual Real operator()() = 0; // function value at current x 00039 // set current x 00040 virtual void Set(Real X); // set x, check OK 00041 Real operator()(Real X) { Set(X); return operator()(); } 00042 // set x, return value 00043 virtual bool IsValid(Real X); 00044 operator Real(); // implicit conversion 00045 virtual ~R1_R1() {} // keep gnu happy 00046 }; 00047 00048 class SolutionException : public BaseException 00049 { 00050 public: 00051 static unsigned long Select; 00052 SolutionException(const char* a_what = 0); 00053 }; 00054 00055 class OneDimSolve 00056 { 00057 R1_R1& function; // reference to the function 00058 Real accX; // accuracy in X direction 00059 Real accY; // accuracy in Y direction 00060 int lim; // maximum number of iterations 00061 00062 public: 00063 OneDimSolve(R1_R1& f, Real AccY = 0.0001, Real AccX = 0.0) 00064 : function(f), accX(AccX), accY(AccY) {} 00065 // f is an R1_R1 function 00066 Real Solve(Real Y, Real X, Real Dev, int Lim=100); 00067 // Solve for x in Y=f(x) 00068 // X is the initial trial value of x 00069 // X+Dev is the second trial value 00070 // program returns a value of x such that 00071 // |Y-f(x)| <= accY or |f.inv(Y)-x| <= accX 00072 00073 private: 00074 Real x[3], y[3]; // Trial values of X and Y 00075 int L,C,U,Last; // Locations of trial values 00076 int vpol, hpol; // polarities 00077 Real YY; // target value 00078 int i; 00079 void LookAt(int); // get new value of function 00080 bool Finish; // true if LookAt finds conv. 00081 bool Captured; // true when target surrounded 00082 void VFlip(); 00083 void HFlip(); 00084 void Flip(); 00085 void State(int I, int J, int K); 00086 void Linear(int, int, int); 00087 void Quadratic(int, int, int); 00088 }; 00089 00090 00091 #ifdef use_namespace 00092 } 00093 #endif 00094 00095 // body file: solution.cpp 00096 00097 00099