c - 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 被设计为迭代调用以接近正确答案。您不能直接使用 20 作为参数 n 来调用它。您必须先用 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 中。也许这有一些教学上的借口，但这样的代码应该永远在现实世界的代码中使用。 (实现进行迭代优化的例程的正确方法是将状态信息返回给调用者，他们可以在以后的调用中传回这些信息。)

调用函数的说明在您引用的引用资料中:

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.

不要在没有理解的情况下盲目复制源码。