前言

前言

simpson积分

simpson积分公式
$${\int }_a^{b}f(x)dx\approx \frac{b-a}{6}[f(a)+f(b)+4f(\frac{a+b}{2})]$$
与梯形积分类似，当区间[a,b]较大时，积分与实际相差较大。

复化simpson积分

假设区间[a,b]被划分为n个小区间，则有区间间隔为$h=\frac{b-a}{n},x_k=a+kh(k=0,\cdots ,n)$，对于每一个小区间都有：
$${\int }_{x_k}^{x_{k+1}}f(x)dx\approx \frac{h}{6}[f(x_k)+f(x_{k+1})+4f(x_{k+\frac{1}{2}})]$$
在[a,b]区间内进行累加可以得到
$${\int }_a^{b}f(x)dx\approx \frac{h}{6}[f(a)+f(b)+2\sum _{k=0}^{n-1}f(x_{k+1})+4\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}})]=S_n$$

误差分析

复化simpson积分的余项为：$$\begin{gathered} R[f]=-\frac{b-a}{180}(\frac{h}{2}{)}^4{f}^{(4)}(\xi ) \\ =-\frac{(b-a{)}^5}{2880n^4}{f}^{(4)}(\xi ) \end{gathered}$$
与梯形积分相似，对其加密一倍n，可以得到对应的事后误差估计公式：
$$S_{2n}(f)-s_n(f)<15ϵ$$
可以等价于：
$$I(f)-S_{2n}(f)<ϵ$$

代码

#include<iostream>
#include<cmath>
#include<iomanip>//这个头文件仅仅是用来设置cout的输出精度
#define abs(x) (x>0?x:-x)
using namespace std;
double Simpson(int n,double a,double b,float(*f)(float x))
{
    double f_x=0.0f;
    double h=(b-a)/n;

    for(int i=0;i<n;i++)
    {
        f_x+=4*f(a+i*h+h/2);//算f(x_(k+1/2))
    }
    for(int i=1;i<n;i++)
    {
        f_x+=2*f(a+i*h);//算f(x_(k+1))
    }
    f_x+=f(a)+f(b);
    f_x=f_x*h/6;
    return f_x;
}
//直接在这里换函数，函数为sin(x)/x
float fx(float x)
{
    float result;
    float x_temp=((x==0)?1e-15:x);
    result=sin(x_temp)/x_temp;
    return result;
}
int main()
{
    double error=1e-6;//表示误差小于10^-6次方
    double a = 0.0, b = 2.0;
    int n=1;
    double f_x_n=(fx(a)+fx(b))*(b-a)/2;
    double f_x_2n;
    while(true)
    {
        f_x_2n=Simpson(n*2,a,b,fx);//算T2n
        if(fabs(f_x_n-f_x_2n)<(error*15))
        {
            // cout<<n<<":simposon error="<<fabs(f_x_n-f_x_2n)<<endl;
            // printf("%.8f,%.8f\n",f_x_n,f_x_2n);
            cout << "Simpson积分的误差为：" << fabs(f_x_n - f_x_2n) << endl;
            n*=2;
            break;
        }
        n+=1;
        f_x_n=Simpson(n,a,b,fx);
    }
    cout<<"区间划分数量:n="<<n<<"，积分值为："<<std::setprecision(10)<<f_x_2n<<endl;
    return 0;
}

结果分析

与梯形积分相比，在相同的被积函数和积分区间上，为了达到一样的误差范围，其迭代次数更少：

使用matlab进行对比可以得到：

可以看出其实际误差已经达到了期望。