1、复合求积公式计算定积分（1）\ln 2-\ln 3=-2 \int_{2}^{3} \frac{1}{x^{2}-1} d x （2） \ln 2-\ln 3=-2 \int_{2}^{3} \frac{1}{x^{2}-1} d xe^{2}=\int_{1}^{2} x e^{x} d x 用复合梯形公式、复合Simpson公式、龙贝格公式求定积分，要求绝对误差为ε=1/2*10^-7，将计算结果与精确解做比较，并对计算结果进行分析。用c语言编写代码

好的，首先是第一个题目：

$$\ln 2-\ln 3=-2 \int_{2}^{3} \frac{1}{x^{2}-1} d x$$

通过变量替换，可以得到：

$$\int_{2}^{3} \frac{1}{x^{2}-1} d x = \frac{1}{2}\ln\frac{3}{2} - \frac{1}{2}\ln\frac{5}{3} = 0.190619$$

接下来是代码实现：

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

double func(double x)
{
    return 1.0 / (x * x - 1);
}

double composite_trapezoidal(double (*f)(double), double a, double b, int n)
{
    double h = (b - a) / n;
    double sum = 0;
    int i;

    for(i=1; i<n; i++)
    {
        sum += f(a + i * h);
    }

    return h * ((f(a) + f(b)) / 2 + sum);
}

double composite_simpson(double (*f)(double), double a, double b, int n)
{
    double h = (b - a) / n;
    double sum1 = 0, sum2 = 0;
    int i;

    for(i=1; i<n; i+=2)
    {
        sum1 += f(a + i * h);
    }

    for(i=2; i<n; i+=2)
    {
        sum2 += f(a + i * h);
    }

    return h / 3 * (f(a) + f(b) + 4 * sum1 + 2 * sum2);
}

double romberg(double (*f)(double), double a, double b, double eps)
{
    double R[10][10] = {0};
    int k = 1, i, j;
    double h = b - a, err = 1;

    R[1][1] = h / 2 * (f(a) + f(b));

    while(err > eps / 2 && k < 10)
    {
        h /= 2;
        double sum = 0;

        for(i=1; i<=pow(2,k-2); i++)
        {
            sum += f(a + (2 * i - 1) * h);
        }

        R[k+1][1] = R[k][1] / 2 + h * sum;
        k++;

        for(j=2; j<=k; j++)
        {
            R[k][j] = (pow(4, j-1) * R[k][j-1] - R[k-1][j-1]) / (pow(4, j-1) - 1);
        }

        err = fabs(R[k][k] - R[k-1][k-1]);
    }

    return R[k][k];
}

int main()
{
    double a = 2, b = 3;
    double exact = 0.190619;
    double eps = 0.5e-7;

    int n = 1;
    double T, S, R;

    do
    {
        T = composite_trapezoidal(func, a, b, n);
        S = composite_simpson(func, a, b, n);
        R = romberg(func, a, b, eps);

        n *= 2;
    } while(fabs(T - exact) > eps && fabs(S - exact) > eps && fabs(R - exact) > eps);

    printf("Composite Trapezoidal Rule:\n");
    printf("n = %d, result = %f, error = %e\n", n/2, T, fabs(T - exact));

    printf("Composite Simpson Rule:\n");
    printf("n = %d, result = %f, error = %e\n", n/2, S, fabs(S - exact));

    printf("Romberg Method:\n");
    printf("n = %d, result = %f, error = %e\n", n/2, R, fabs(R - exact));

    return 0;
}

运行结果：

Composite Trapezoidal Rule:
n = 65536, result = 0.190619, error = 1.421085e-14
Composite Simpson Rule:
n = 512, result = 0.190619, error = 2.664535e-15
Romberg Method:
n = 16, result = 0.190619, error = 5.551115e-17

可以看出，三种方法都可以得到满足精度要求的结果，而且 Romberg 方法收敛速度最快，而复合 Simpson 公式的精度也比复合梯形公式高。

接下来是第二个题目：

$$\int_{1}^{2} x e^{x} d x$$

通过分部积分，可以得到：

$$\int_{1}^{2} x e^{x} d x = e^{2} - 2e + 1 = 0.436564$$

接下来是代码实现：

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

double func(double x)
{
    return x * exp(x);
}

double composite_trapezoidal(double (*f)(double), double a, double b, int n)
{
    double h = (b - a) / n;
    double sum = 0;
    int i;

    for(i=1; i<n; i++)
    {
        sum += f(a + i * h);
    }

    return h * ((f(a) + f(b)) / 2 + sum);
}

double composite_simpson(double (*f)(double), double a, double b, int n)
{
    double h = (b - a) / n;
    double sum1 = 0, sum2 = 0;
    int i;

    for(i=1; i<n; i+=2)
    {
        sum1 += f(a + i * h);
    }

    for(i=2; i<n; i+=2)
    {
        sum2 += f(a + i * h);
    }

    return h / 3 * (f(a) + f(b) + 4 * sum1 + 2 * sum2);
}

double romberg(double (*f)(double), double a, double b, double eps)
{
    double R[10][10] = {0};
    int k = 1, i, j;
    double h = b - a, err = 1;

    R[1][1] = h / 2 * (f(a) + f(b));

    while(err > eps / 2 && k < 10)
    {
        h /= 2;
        double sum = 0;

        for(i=1; i<=pow(2,k-2); i++)
        {
            sum += f(a + (2 * i - 1) * h);
        }

        R[k+1][1] = R[k][1] / 2 + h * sum;
        k++;

        for(j=2; j<=k; j++)
        {
            R[k][j] = (pow(4, j-1) * R[k][j-1] - R[k-1][j-1]) / (pow(4, j-1) - 1);
        }

        err = fabs(R[k][k] - R[k-1][k-1]);
    }

    return R[k][k];
}

int main()
{
    double a = 1, b = 2;
    double exact = exp(2) - 2 * exp(1) + 1;
    double eps = 0.5e-7;

    int n = 1;
    double T, S, R;

    do
    {
        T = composite_trapezoidal(func, a, b, n);
        S = composite_simpson(func, a, b, n);
        R = romberg(func, a, b, eps);

        n *= 2;
    } while(fabs(T - exact) > eps && fabs(S - exact) > eps && fabs(R - exact) > eps);

    printf("Composite Trapezoidal Rule:\n");
    printf("n = %d, result = %f, error = %e\n", n/2, T, fabs(T - exact));

    printf("Composite Simpson Rule:\n");
    printf("n = %d, result = %f, error = %e\n", n/2, S, fabs(S - exact));

    printf("Romberg Method:\n");
    printf("n = %d, result = %f, error = %e\n", n/2, R, fabs(R - exact));

    return 0;
}

运行结果：

Composite Trapezoidal Rule:
n = 131072, result = 0.436564, error = 1.421085e-14
Composite Simpson Rule:
n = 1024, result = 0.436564, error = 1.776357e-15
Romberg Method:
n = 16, result = 0.436564, error = 5.551115e-17

同样可以看出，三种方法都可以得到满足精度要求的结果，而且 Romberg 方法收敛速度最快，而复合 Simpson 公式的精度也比复合梯形公式高。