home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
OS/2 Shareware BBS: 10 Tools
/
10-Tools.zip
/
autopsp.zip
/
REGRESS
/
CSIMP.CPP
< prev
next >
Wrap
C/C++ Source or Header
|
1995-07-22
|
3KB
|
113 lines
// CSimpson
// Method definition for CSimpson integration class
#include "csimp.h"
CSimpson::CSimpson(void)
{
n=20; // Safe values.
result=0; // Ensure no uninitialization problems
cur_result=0;
old_result=0;
result=0;
width=0;
tolerance=0.005; // Resonable tolerance.
}
// I couldn't compile with out this place holder.
double CSimpson::IntMe(double x)
{
return(x);
}
// This method handles the infinities and can add on the
// "other half" if the user passes in a non-zero total.
// Total should be the total area under the curve if known.
double CSimpson::Integrate(double low, double high, double total)
{
result=0;
if (low==INFINITY && high==INFINITY){
result=total; // Easy case first
return(result);
}
else if (low==INFINITY && high<=0){
CalcSimpsons(high,0);
result += (total/2) - cur_result;
}
else if (low==INFINITY && high>0){
CalcSimpsons(0,high);
result += (total/2) + cur_result;
}
else if (low<=0 && high==INFINITY){
CalcSimpsons(low,0);
result += (total/2) + cur_result;
}
else if (low>0 && high==INFINITY){
CalcSimpsons(0,low);
result += (total/2) - cur_result;
}
else {
CalcSimpsons(low,high);
result=cur_result;
}
return(result);
}
// Calculates the integral using simpsons rule.
void CSimpson::CalcSimpsons(double low, double high)
{
double sum=0;
double slider; // Slides across the function.
if(low==high){ // Handles a special case.
cur_result=0; // Space between x and x is 0.
return;
}
old_result=INFINITY; // Avoid accidental immediate kick out.
n=10; // Starting at 10 makes inital run have 20 blocks.
do{
sum=0; // Summation starts at 0 for each attempt.
n*=2; // double number of blocks.
old_result=cur_result; // Update old_result before changing cur_result.
width=(high-low)/n; // Calc the width.
slider=0; // Slide starts at 0.
sum+=IntMe(low+slider); // Following loop implements simpsons algorythm.
slider+=width;
sum+=4*IntMe(low+slider);
while(low+slider <= high-(2*width)){
slider+=width;
sum+=2*IntMe(low+slider);
slider+=width;
sum+=4*IntMe(low+slider);
}
sum+=IntMe(high);
cur_result=(width/3)*sum; // Simpsons fudge factor.
// Loop until the value is sufficiently
// cool to be accepted into the class.
} while(fabs(cur_result-old_result) >= tolerance);
}
// Calculates the Gamma function based on the books description.
// Will only handle integers and (integers + .5)...
// (Accepts the set of all reals divisible by .5?)
double CSimpson::Gamma(double x)
{
if (x>1){
return((x-1)*Gamma(x-1));
}
else if (x==.5){
return(pow(M_PI,.5));
}
return(1);
}
// Sets the private member: tolerance.
void CSimpson::SetTolerance(double Tol)
{
tolerance=Tol;
}