Clever Geek Handbook

Boole's rule for N intervals (C)

I am trying to implement a Boole rule over n intervals using this formula boole's rulex =

So far I have developed this code:

//f = function on the range [a,b] n = number of intervals
long double booles(long double (*f) (long double), 
                double a, double b, int n) 
{
  long double sum=7*f(a); //because the starting value is not doubled 
  long double h = (b - a)/(n-1); //width of each interval
  int mod;
  int i =1;
  while (i<n-1)
    {
      mod = i%4;
      if (mod==0)
        sum+=14*f(a+i*h);
      else if (mod==1)
        sum+=32*f(a+i*h);
      else if (mod==2)
        sum+=12*f(a+i*h);
      else
        sum+=32*f(a+i*h);
      ++i;
    }
  return 2* h/45 * sum;
}

This will run and give a decent answer, but not in the rules for the Bool error, which indicates an error error. I fixed the problem with the first term doubling, but I'm not sure how to fix the possible doubling at the end of the loop. Also, the error is relatively large and I don't think my only problem is the last four terms.

+3


source to share


1 answer


  • long double

    The wiki says: extended precision floating point type. Unspecified actual properties. Unlike float and double types, it can be either 80-bit floating point or IEEE 754 four-point floating point, but not IEEE 754 if a higher precision format is provided, otherwise it is the same as double. See the long twin article for details.

    • so it's hard to tell which datatype you are using (I prefer to use doubles)

  • <strong> constants

    you are mixing integers and floats together, so the compiler has to decide which to use. Rewrite all floats as 45

    in 45.0

    to make sure it's done correctly, or a+i*h

    ... i

    - int and h

    - double ...

  • integration

    don't know the magnitude of the sum and the range of your values, but to improve the floating precision, you should avoid adding large and small values ​​together , because if the exponents are too different, you lose too many corresponding bits of the mantissa.

    So the sum in two variables is something like this (in C ++):

    double suml=0.0,sumh=0.0,summ=1000.0;
for (int i=0;i<n;i++)
 {
 suml+=...; // here goes the formula
 if (fabs(suml)>summ) { sumh+=suml; suml=0; }
 } sumh+=suml; suml=0;
// here the result is in sumh

    • summ

      is the maximum value of suml. It should be in a relatively safe range compared to your sum iteration value, for example 100-10000

      times the average.

    • suml

      is the variable summation of small quantities

    • sumh

      - variable sum of large values

    If the range of your summed values ​​is really large, you can add another if

    if (fabs(value)>summ) sumh+=value; else suml+=value;

    if it is even much larger than you can sum to any number of variables in the same way, just divide the range of values ​​by some sense of the full ranges

  • Formula

    Maybe I'm missing something, but why are you moderating? As I can see you don't need a loop at all and also ifs are deprecated, so why use a+i*h

    instead a+=h

    ? this will improve performance and accuracy.

    something like that:

    double sum,h;
h = (b-a)/double(n-1);
sum = 7.0*f(a); a+=h;
sum+=32.0*f(a); a+=h;
sum+=12.0*f(a); a+=h;
sum+=32.0*f(a); a+=h;
sum+= 7.0*f(a); a+=h;
return 2.0*h*sum/45.0;
// + the thing in the bullet 3 of coarse ...
// now I see you had an error in your constants !!!

[edit1] Done (not quadrupled)



//---------------------------------------------------------------------------
double f(double x)
    {
//  return x+0.2*x*x-0.001*x*x*x+2.0*cos(0.1*x)*tan(0.01*x*x)+25.0;
    return x+0.2*x*x-0.001*x*x*x;
    }
//---------------------------------------------------------------------------
double integrate_rect(double (*f)(double),double x0,double x1,int n)
    {
    int i;
    double s=0.0,x=x0,dx=(x1-x0)/double(n-1);
    for (i=0;i<n;i++,x+=dx) s+=f(x);
    return s*dx;
    }
//---------------------------------------------------------------------------
double integrate_Boole(double (*f)(double),double x0,double x1,int n)
    {
    n-=n%5; if (n<5) n=5;
    int i;
    double s=0.0,x=x0,dx=(x1-x0)/double(n-1);
    for (i=0;i<n;i+=5)
        {
        s+= 7.0*f(x); x+=dx;
        s+=32.0*f(x); x+=dx;
        s+=12.0*f(x); x+=dx;
        s+=32.0*f(x); x+=dx;
        s+= 7.0*f(x); x+=dx;
        }
    s*=(2.0*dx)/(45.0);
    return s*1.25; // this is the ratio to cover most cases
    }
//---------------------------------------------------------------------------
void main()
    {
    int n=100000;
    double x0=0.0,x1=+100.0,i0,i1;
    i0=integrate_rect (f,x0,x1,n); cout << i0 << endl;
    i1=integrate_Boole(f,x0,x1,n); cout << i1 << endl << i0/i1;
    }
//---------------------------------------------------------------------------

I use the rectangle rule because FPU is the fastest and most accurate way. More advanced approaches are better on paper, but on a computer, the added overhead and rounding is usually slower / less accurate than the same precision with a rectangular rule.

+5


source






More articles:


All Articles