使用 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 的值不同...我期待

编辑：我现在对上面的代码进行了建议的更改，这要感谢下面几个用户的 cmets。非常感谢您的帮助！我现在的输出有问题，请参阅终端的附加图像。对于此输入，它应该在 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)，而对于这种类型的积分，您应该考虑增加上限，并且可能是不均匀的网格间距。