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

Question

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

#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 中。也许有一些教学上的借口，但这样的代码应该永远在现实世界的代码中使用。 （实现进行迭代细化的例程的适当方法是将状态信息返回给调用者，它们可以在未来的调用中传回。）

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

当以 n=1 调用时，例程返回来自 a 的 [f(x) dx 的积分的最粗略估计到b]。后续调用 n=2,3,... a（按该顺序）将通过添加 2^n-2 个额外的内部点来提高准确性。

不要在不懂的情况下盲目复制源代码。