C语言实现Trapezoid Rule（梯形法则）输出错误，求排查

Let's fix your code and break down the issues clearly—you're just one critical mistake away from getting the right results.

The Root Cause: Hardcoded Function Instead of Using the Function Pointer #

Your trapzd function accepts a function pointer func as its first parameter, but you never actually use it! Everywhere in the function, you're calling h(x) directly, regardless of what function you pass in from main. This means:

• When you try to integrate sin(x), it technically uses the right function, but there's a subtle initialization quirk that throws off early results.
• When you switch to exp(x), the code still calculates the integral of sin(x) instead of exp(x)—hence the completely wrong output.

Fixing the Code #

Replace every instance of h(a), h(b), and h(x) in trapzd with func(a), func(b), and func(x) respectively. We've also added a wrapper for exp to match the function signature expected by trapzd:

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

float h(float x) { return sin(x); }
float exp_func(float x) { return exp(x); } // Wrap standard exp to match float signature

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) * (func(a) + func(b))); // Use passed-in func instead of hardcoded h
    } 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 += func(x); // Use func here too
        s = 0.5 * (s + (b - a) * sum / tnm);
        return s;
    }
}

int main(void) {
    // Test sin(x) from 0 to 1 (expected ~0.4597)
    printf("Integrating sin(x) from 0 to 1:\n");
    for (int i = 1; i <= 10; ++i) {
        printf("With n = %d, the approximation is %g.\n", i, trapzd(h, 0, 1, i));
    }

    // Test exp(x) from 0 to 1 (expected ~1.7183)
    printf("\nIntegrating exp(x) from 0 to 1:\n");
    for (int i = 1; i <= 10; ++i) {
        printf("With n = %d, the approximation is %g.\n", i, trapzd(exp_func, 0, 1, i));
    }
}

Additional Notes #

1. Function Signature Matching: The standard exp function takes and returns double, so we wrapped it in exp_func to match the float(*func)(float) signature expected by trapzd. For better precision in numerical calculations, you could also modify trapzd to use double throughout.
2. Static Variable s: The static variable is intentional—it preserves the previous approximation to avoid redundant calculations, which is part of the efficient trapezoid rule implementation from Numerical Recipes. Just note that if you call trapzd for different integrals without resetting s, it will carry over old values. You could add a reset mechanism or pass s by reference if needed.

Expected Results #

After fixing, you'll see the approximation for sin(x) converge to ~0.4597, and exp(x) converge to ~1.7183 as n increases—matching your expected values perfectly.

内容的提问来源于stack exchange，提问作者Eren