【数值分析】复化积分公式

对于积分$\int_{a}^{b}f(x)dx$只要找到被积公式的原函数$F(x)$，利用牛顿莱普利兹公式有：

$$I=\int_{a}^{b}f(x)dx=F(b)-F(a)$$

但是，实际使用这种求积分的方法往往是有困难的，因为大量的被积函数的原函数是不能用初等函数表示的；另外，当$f(x)$是由测量或数值计算给出的一张数据表时，牛顿莱普利兹公式也无法直接运用，因此有必要研究积分的数值计算问题。

---

对于一些理论的推导，大家可以看看维基百科，下面我主要给出牛顿-科特斯公式在$n=1$(梯形求积公式)、$n=2$（辛普森公式）的情况，并通过代码实现。

梯形公式：

$$T=(b-a)\left[\frac{1}{2}f(a)+\frac{1}{2}f(b)\right]$$

辛普森公式：

$$S=(b-a)\left[\frac{1}{6}f(a)+\frac{4}{6}f(\frac{a+b}{2})+\frac{1}{6}f(b)\right]$$

应用高阶牛顿-科特斯公式计算积分时，会出现数值不稳定的情况，而低阶公式往往因为积分步长过大使得离散误差变大，因此，为了提高求积公式的精度，可以把积分区间分成若干个子区间，在每个子区间上使用低阶求积公式，然后将结果加起来，这种方法称为复化求积法。

复化梯形公式

将区间$[a,b]$划分为$n$等分，步长为$h=(b-a)/h$，节点为$xi=a+ih,i=1,2,\cdots,n+1$，在每个子区间$[x_i,x{i+1}]$ 使用梯形公式得：

$$\int_{a}^{b}f(x)dx\approx\sum\limits_{i=1}^{n}\int_{x_i}^{x_{i+1}}f(x)dx=T_n=\frac{h}{2}\left[f(a)+2\sum\limits_{i=1}^{n-1}f(a+ih)+f(b)\right],h=\frac{b-a}{n}$$

复化辛普森公式

根据复化梯形公式的推导，同理可得复化辛普森公式为：

$$S_n=\frac{h}{6}\sum\limits_{i=0}^{n-1}\left[f(x_i)+4f(x_{i+1/2})+f(x_{i+1})\right], h=\frac{b-a}{n}$$

---

下面我们通过实例来实现复化梯形公式和复化辛普森公式：

对于函数$f(x)=\sin(x)/x$，试用复化梯形公式和复化辛普森公式计算函数$f(x)$在$[0, 1]$上的积分。

具体的程序实现如下：

#include<stdio.h>
#include<math.h>
double Function(double x)//所要计算积分的函数f(x)
{
    if(x==0)//sin(x)/x在0处的取值为1
        return 1;
    else
        return sin(x)/x;
}
//复化梯形公式
double Trapz(double a,double b,int n)
{
    double h=(b-a)/n;
    double T=0;
    for(int i=1;i<n;i++)
    {
        T=T+Function(a+i*h);
    }
    T*=2;
    T=(Function(a)+Function(b)+T)*h/2;
    return T;
}
//复化辛普森公式
double MulripleSimpson(double a,double b,int n)
{
    double h=(b-a)/n;
    double T=0;
    for(int i=0;i<n;i++)
    {
        T=T+Function(a+i*h)+4*Function(a+(i+0.5)*h)+Function(a+(i+1)*h);
    }
    T=T*h/6;
    return T;
}
void main()
{
    printf("使用复化梯形公式可得：%f\n",Trapz(0,1,8));
    printf("使用复化辛普森公式可得：%f\n",MulripleSimpson(0,1,4));
}

运行结果如下图：

结果分析：

比较复化梯形公式和复化辛普森公式两种方法的运行结果，我们发现复化辛普森公式与准确值$0.9460831$更加接近，复化梯形公式只有$2$位有效数字，而复化辛普森公式有$6$为有效数字。