c - c集成程序中的舍入错误

Question

我正在编写一个程序，针对不同的幂 n 值和常数 a 计算函数的积分。我的程序似乎运行正常，但我的结果出现了一个小的舍入错误，我不知道为什么。我知道我有一个错误，因为我的一个 friend 也在制作相同的程序，他的结果与我的略有不同，他的结果绝对是正确的，因为在计算器上进行积分给出的值更接近他的。下面是我的结果和他对 a=2 和 n=1 的结果。

他的结果:0.189070
我的结果:0.189053

我尝试过遍历和类型转换我能想到的几乎所有东西，但仍然无法弄清楚我从哪里得到我的错误，任何帮助指出我是一个白痴的地方将不胜感激! :p

我的程序:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define debug 0
#define N (double)10000

double Integrand(double x, int a, int n);
double Integral(double *x, double dx, int a, int n);

int main (int argc, char* argv[])
{
    int j,a,n=0,count=0,size=(int)N;
    double dx=1/N, x[size];

    sscanf(argv[1], "%d", &a);
    for(j=0;j<N;j++) {
        x[j]=(double)(j)*dx;
    }
    for(n=1;n<=10;n++) {
        printf("n is %d integral is %lf\n",n,Integral(x,dx,a,n));
    }    
    return(EXIT_SUCCESS);
}

double Integral(double *x, double dx, int a, int n)
{
    int i;
    double result=0;

    for(i=0;i<N;i++) {
       result +=(double)((Integrand((double)x[i],a,n))*dx);
    }
    return(result);
}

double Integrand(double x, int a, int n)
{
   double result;
   result=(double)(((pow(x,(double)n))/(x+(double)a)));
   return(result);
}

Answer 1

这不是舍入错误，您只是没有为积分点选择最佳选择。将初始化更改为

x[j]=(j+0.5)*dx;

这样您就可以取每个积分条的中点来计算被积函数的值。如果您总是取左端点或右端点，对于单调函数，您将得到系统性的过大误差。

如果你用黎曼和来近似一个足够平滑的函数 f 的积​​分，

 b           n
 ∫ f(x) dx ≈ ∑ f(y_k)*(b-a)/n
 a          k=1

区间[x_(k-1), x_k] = [a+(k-1)*(b-a)/n, a+k*( b-a)/n]影响误差和收敛速度。写作

f(x) = f(y_k) + f'(y_k)*(x-y_k) + 1/2*f''(y_k)*(x-y_k)² + O((x-y_k)³)

在那段时间里，你发现

x_k                                   x_k                            x_k
 ∫ f(x) dx = f(y_k)*(b-a)/n + f'(y_k)* ∫ (x-y_k) dx + 1/2*f''(y_k) * ∫ (x-y_k)² dx + O(1/n^4)
x_(k-1)                               x_(k-1)                       x_(k-1)

           = f(y_k)*(b-a)/n + 1/2*f'(y_k)*(b-a)/n*((x_k-y_k)-(y_k-x_(k-1))) + O(1/n³)

关于近似值 f(y_k)*(b-a)/n 的第一个和最大误差项消失了

y_k = (x_k + x_(k-1))/2

为您提供该 strip 的总误差 O(1/n³)，以及整个黎曼和的总误差 O(1/n²)。

如果您选择y_k = x_(k-1)(或y_k = x_k)，第一个误差项变为

±1/2*f'(y_k)*[(b-a)/n]²

导致 O(1/n) 的总误差。