0

我正在制作一个程序，该程序针对不同的幂 n 值和常数 a 找到函数的积分。我的程序似乎工作正常，但我的结果中有一个小的舍入错误，我不知道为什么。我知道我有一个错误，因为我的一个朋友也在制作相同的程序，他的结果与我的略有不同，而且他绝对是正确的，因为在计算器上进行积分可以得到更接近他的值。下面是我的结果和他的 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);
}
4

2 回答 2

3

It's not a rounding error, you just don't pick the best choice for integration points. Change the initialisation to

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

so that you take the midpoint of each integration strip to calculate the value of the integrand. If you always take the left or right endpoint, you will get a systematically too large error for monotonic functions.

If you approximate the integral of a sufficiently smooth function f by a Riemann sum,

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

the choice of y_k in the interval [x_(k-1), x_k] = [a+(k-1)*(b-a)/n, a+k*(b-a)/n] influences the error and the speed of convergence. Writing

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

in that interval, you find that

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³)

and the first and largest error term with respect to the approximation f(y_k)*(b-a)/n vanishes for

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

giving you an overall O(1/n³) error for that strip, and a total O(1/n²) error for the entire Riemann sum.

If you choose y_k = x_(k-1) (or y_k = x_k), the first error term becomes

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

leading to an O(1/n) total error.

于 2012-12-01T11:03:53.133 回答
0

在 Linux 终端提示符下，键入：

man fegetenv
于 2012-12-01T10:59:06.337 回答