C中的梯形规则，错误答案

Question

我想用C语言为一个项目实现梯形元组。我正在使用第 137 页 Numerical Recipes in C 一书中的示例。我首先只想复制代码，尝试不同的功能，然后尝试理解每一行。但是，我的程序没有输出正确的答案：

#include <stdio.h>
#include <math.h>
float h(float x)
{ 
    return sin(x);
}
float trapzd(float (*func)(float), float a, float b, int n)
{   float x, tnm, sum, del; 
    static float s;
    int it, j;
    if (n == 1) {
        return (s=0.5*(b-a)*(h(a)+h(b)));


    } 
    else
    {
        for (it=1,j=1;j<n-1;j++) it <<=1;
        tnm=it;
        del=(b-a)/tnm;
        x=a+0.5*del;
        for(sum=0.0,j=1; j<=it; j++, x+=del) sum+=h(x);
        s=0.5*(s+(b-a)*sum/tnm);   

        return s;
    }

}
int main(void)
{   
    for (int i = 1; i <= 10; ++i)
    {
    printf("With n = %d, the approximation is %g.\n", i, trapzd(h,0,1,i));
    }

}

所以我想将 sin(x) 从 0 整合到 1，这应该等于 0.45 但程序输出 0.000002。如果我使用 exp(x)，输出是 0.50000 而不是大约 1.7。为什么这不起作用？

Answer 1

例程 trapzd 旨在迭代调用以接近正确答案。您不能直接使用参数 n 的 20 调用它。您必须先用 1 调用它，然后是 2，然后是 3，依此类推。

例如，使用 Numerical Recipes 中的原始代码（删除您在 trapzd 中插入的 printf 语句）：

for (int i = 1; i <= 20; ++i)
    printf("With n = %d, the approximation is %g.\n", i, trapzd(h, 0, 1, i));

此外，正如 Mark Dickinson 指出的那样，x=+del 应该是 x+=del。

Numerical Recipes 中的来源是可怕的 设计。该函数具有内部状态，在静态对象 s 中。也许这有一些教学上的借口，但这样的代码应该 永远不会 在现实世界的代码中使用。（实现进行迭代优化的例程的正确方法是 return 向调用者声明信息，他们可以在以后的调用中传回这些信息。）

调用函数的说明在您引用的参考文献中：

When called with n=1, the routine returns the crudest estimate of [the integral of f(x) dx from a to b]. Subsequent calls with n=2,3,... a (in that sequential order) will improve the accuracy by adding 2^n-2 additional interior points.

不要在没有看懂的情况下盲目复制源码

C中的梯形规则，错误答案

Trapezoid Rule in C, wrong answer

c

numerical-integration