c++ - N 个区间的 bool 规则 (C)

Question

我正在尝试使用此公式 在 n 个区间上实现 bool 规则

到目前为止，我已经开发了这段代码:

//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;
}

这将运行并给出一个体面的答案，但它不在 Bool 错误的规则之内，该规则指出错误是 。 .我解决了将第一项加倍的问题，但我不确定如何解决循环结束时可能加倍的问题。此外，误差相对较大，我不认为我唯一的问题是最后四个术语。

Answer 1

1. 长双

   Wiki 说: 扩展精度浮点类型。实际属性未指定。与 float 和 double 类型不同，它可以是 80 位浮点格式，非 IEEE“double-double”或 IEEE 754 四精度浮点格式(如果提供了更高精度的格式)，否则与双倍的。有关详细信息，请参阅有关 long double 的文章。

   • 所以很难说你实际使用的是什么数据类型(我更喜欢使用 double )

2. 常量

   您将整数和 float 混合在一起，因此编译器必须决定使用什么。将所有 float 重写为 45 到 45.0 以确保它正确完成或 a+i*h ... i 是 int 而 h 是 double ...

3. 整合

   不知道值的总和和范围的大小，但要提高 float 精度，您应该避免将大值和小值相加，因为如果指数差异太大，您就会失去太多相关信息尾数位。

   所以两个变量的总和是这样的(在 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 是 suml 的最大值。与迭代总和值相比，它应该处于相对安全的范围内，例如 100-10000 是平均值的倍数。

   • suml 是低量级和变量

   • sumh 是大数量级的和变量

   如果你求和的范围真的很大那么你可以添加另一个如果

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

   如果它更大，那么您可以以相同的方式对任意数量的变量求和，只需将值的范围划分为一些有意义的完整范围即可

4. 公式

   可能是我遗漏了什么，但你为什么要修改？如我所见，您根本不需要循环，而且 ifs 也已过时，那么为什么要使用 a+i*h 而不是 a+=h？它将提高性能和精度

   像这样:

   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] 部门实现(不是四倍)

//---------------------------------------------------------------------------
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;
    }
//---------------------------------------------------------------------------

我主要使用矩形规则，因为在 FPU 上是最快和最精确的方法。更高级的方法在纸面上更好，但在计算机上，增加的开销和舍入通常比矩形规则的相同精度更慢/更不精确