c集成程序中的舍入错误

Posted 2023-02-25

【中文标题】c集成程序中的舍入错误

【英文标题】：Rounding error in c integration program

【发布时间】：2012-12-01 10:53:59

【问题描述】：

我正在编写一个程序，该程序可以针对不同的幂 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);

【问题讨论】：

你们使用的是同一个编译平台吗？你读过***.com/q/13571073/139746 吗？

您的好友代码与您自己的代码相同吗？

而不是直接使用 dx = 1/N 计算 1/N。我认为当你在 dx 中取 1/N 的值时，会有一些精度损失。它正在传播。

“我尝试过并铸造几乎所有我能想到的东西”……你在那里学到了宝贵的一课：只铸造你能想到的所有东西并不能解决问题，它使程序看起来很可怕。

【参考方案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]中选择y_k会影响误差和收敛速度。写作

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

为您提供该条带的总体 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) 总误差。

【讨论】：

您是否有工具可以创建您粘贴的格式化方程式或者是手工制作的？

@JonasWielicki 所有手工作品。例如，我键入 html 实体名称 ∑，复制显示的字形并将其粘贴到代码块中。 （男孩，我多么希望这个网站支持 LaTeX 渲染。）

感谢您的解释，我现在感觉有点像个白痴，因为我应该意识到，我自己记得几年前在大学里做中点与结束 reimann 求和，现在你已经说过了

@user1831711 如果忘记一个人曾经学过的东西会让一个人成为白痴，我会争夺世界冠军。现在你已经复习了，所以在接下来的几年里，你会记得的。

@user1831711 接受他的回答以感谢他的更新；）

【参考方案2】：

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

man fegetenv

【讨论】：

大多数操作系统都指定新进程以最近舍入的方式启动。我相信 IEEE 754 中甚至有这样的建议。