复化积分的C++实现

分享一下2020年11月11日数值分析实验的解答：

复化梯形公式（Composite Trapezoidal Quadrature）

　　区间[a,b]划分为n等分，分点为$x_k=a+kh， h=\frac{b-a}{n}, k=0,1,2,\dots,n$。对于每个子区间$[x_k, x_{k+1}]$利用公式

$I=\int_{a}^{b}f(x)dx = \sum_{k=0}^{n-1}\int_{x+k}^{x_{k+1}}f(x)dx=\frac{h}{2}\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})]+R_n(f)$

　　复化梯形公式也就是其中的主要部分$T_n=\frac{h}{2}\sum_{k=0}^{n-1}[f(x_k)+f(x_{k+1})]=\frac{h}{2}[f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b)]$

知道这个代码就很容易了

double compositeTrapedoizalQuasdrature(double x[], double y[], int len) {
    double a = x[0];
    double b = x[len - 1];
    double result = function(a) + function(b);
    double h = (b - a) / (len - 1);

    for (int i = 1; i < len-1;i++) {
        result += 2 * function(x[i]);
        printf("i=%d\n", i);
    }

    result *= h * 0.5;
    return result;
}

 

Composite Simpson Quadrature（复化辛普森公式）

　　区间[a,b]划分为n等分，分点为$x_k=a+kh， h=\frac{b-a}{n}, k=0,1,2,\dots,n$。对于每个子区间$[x_k, x_{k+1}]$利用公式

$S=\frac{b-a}{6}f(x)dx=\sum_{k=0}^{n-1}\int_{x_k}^{x_{k+1}}f(x)dz=\frac{h}{6}\sum_{k=0}^{n-1}[f(x_k)+4f(x_{x+\frac{1}{2}})+f(x_{x+1})]+R_n(f)$

复化辛普森公式也就是$S_n=\frac{h}{6}\sum_{k=0}^{n-1}[f(x_k)+4f(x_{k+\frac{1}{2}})+f(x_{k+1})]$

$=\frac{h}{6}[f(a)+4\sum_{k=0}^{n-1}f(x+\frac{1}{2})+2\sum_{k=1}^{n-1}f(x)+f(b)]$

double compositeSimpsonQuasdrature(double x[], double y[], int len) {
    double a = x[0];
    double b = x[len - 1];
    double result = function(a) + function(b);
    double h = (b - a) / (len - 1);

    double xk;
    for (int k = 0; k < len - 1; k++) {
        xk = x[k] + h * 0.5;
        result += 4 * function(xk);
    }

    for (int i = 1; i < len - 1; i++) {
        result += 2 * function(x[i]);
    }

    result *= h;
    result /= 6;
    return result;
}