用 C 编程对 Sin(x)/sqrt(x) 的 0 和无穷大之间的梯形规则进行数值积分

Question

我在 C 编程中找不到 sin(x)/sqrt(x) 的 0 和无穷大之间积分的良好近似值。我正在尝试使用梯形规则。在我的代码中，我还希望用户输入一个精度值，例如 0.001，其中输出值将精确到该小数位数。这段代码出了什么问题？

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

int main(void)
{
int i, N; // integers to interate on in loops //
double h, x, y, precision, lowerLim = 0.0001, upperLim = 10000, f0, fN;

printf("accuracy:");
scanf("%lf", &precision);
N = 10;
double areastored, newarea; 
do {
    newarea = 0.0;
    areastored = newarea; // areastored is the area that I want to compare to the new area calulated as the N (number of partitions) to check the precision of the new value to see if it a better approximation //
    h = (upperLim - lowerLim)/(N-1);
    fN = sin(upperLim)/sqrt(upperLim);
    f0 = sin(lowerLim)/sqrt(lowerLim); // end points evaluated in function //
    newarea = newarea + 0.5*h*(f0 + fN); 

    for (int i = 1; i < N; i++) {
    x = lowerLim + h*i;
    y = sin(x)/sqrt(x);
    newarea = newarea + y*h;   // this loop adds all the middle trapezia areas //
    }

    printf("at N %d integral %f\n", N, newarea);
    N = N*5; // iterate the N so next round of the loop it will approximate an area with more and smaller trapezia //
} while ( fabs ( newarea - areastored ) > precision ); // if this is false then should have an area to the desired precision //

    printf("The integral evaluates to: %lf\n", newarea); 
}

问题是，如果我输入精度 0.01，则计算 N = 10、50、250 的面积，但无法继续，最后一个面积 = 8.53，偏离值 1.253.. . 期待

编辑：由于下面几个用户的评论，我现在已经对上面的代码进行了建议的更改，现在可以看到。非常感谢您的帮助！我现在的输出出现问题，请参阅我的终端的附件图像。对于此输入，它应该在 N 的第 9 次迭代时停止，因此打印的值相当令人费解。为什么会这样？再次感谢您提前提供帮助！

Answer 1

如果你调试你的程序你会看到

h = (upperLim - lowerLim)/N-1;

导致 N=1250 的负值 h。

这会导致死循环

for (x=lowerLim+h; x < upperLim; x+=h) {

因为当 x 达到足够大的绝对值时，当您添加 h 时它将不再改变。

你的意思可能是

h = (upperLim - lowerLim) / (N-1);

Answer 2

您的代码中至少有两个错误

// ...
newarea = 0.0;                       // <-- It's initialized here
do {
    // ...
    h = (upperLim - lowerLim)/N-1;   // <-- This doesn't do what you think
    // ...
    // It's updated, but never reset 
    newarea = newarea + 0.5*h*(f0 + fN);
    // ...
} while ( ... );

您应该将 newarea = 0.0; 行移到循环内并修改计算 h 的公式。

另请注意，您有一个固定的上限 (1000)，而对于此类积分，您应该考虑增加的上限和可能不均匀的网格间距。

用 C 编程对 Sin(x)/sqrt(x) 的 0 和无穷大之间的梯形规则进行数值积分

Numerically integrating with trapezium rule between 0 and infinity of Sin(x)/sqrt(x) with C programming

c

integration

numerical-methods

numerical-integration