复化求积公式

栀梦孤辰

于 2022-04-20 17:36:16 发布

阅读量2.7k

点赞数 7

分类专栏：数值分析文章标签：数据分析算法

本文链接：https://blog.csdn.net/qq_52588437/article/details/124296793

版权

数值分析专栏收录该内容

4 篇文章 1 订阅

订阅专栏

本文介绍了数值积分中提高精度的两种方法：复化梯形公式和复化辛普森公式。通过将区间等分，利用低阶求积公式在每个子区间上计算，再将结果加权求和，从而得到更精确的积分近似值。文章提供了C++代码实现，并以函数f(x)=2/3x^3e^(x^2)为例，对比了两种方法的计算结果和误差。

摘要由CSDN通过智能技术生成

为了提高数值积分的精确度，常将区间 $[a, b]$ 等分成 $n$ 个子区间

其长度为 $h=\frac {b-a}n$ ，
分点为 $x_k = a + kh(k=0,1,...,n)$ ，
在每个子区间上用低阶的求积公式，然后将所有子区间上的计算结果加起来，这样得出的公式称为复化求积公式

这里介绍两个求积公式，复化梯形公式和复化辛普森公式

复化梯形公式

在每个子区间 $x_k,x_{k+1}]$ ，用梯形公式

$\int_{a}^b f(x)dx\approx \frac {b-a}2[f(a)+f(b)]$

得

$T_n = \frac 1n \sum_{k=0}^{n-1} \frac {b-a}2 [ f(x_k) + f(x_{k+1}) ]$

由
$h=\frac {b-a}n \\ x_0=a,x_n=b\\ \sum_{k=0}^{n-1}f(x_k) = f(x_0)+f(x_1)+\cdots+f(x_{n-1}) \\ \sum_{k=0}^{n-1}f(x_{k+1}) = f(x_1)+f(x_2)+\cdots+f(x_n)$

展开化简得

$T_n = \frac h2[ f(a)+2\sum_{k=1}^{n-1}f(x_k)+f(b) ]$

依据上式，使用循环求和，即可实现复化梯形公式
代码如下：

// 复化梯形公式
double CompositeTrapezoidal(double a,double b,double n,double(*f)(double x))
{
    double h = (b - a) / n;    
    double ans = f(a) + f(b);
    for(int i = 1;i < n;i++)
    {
        ans += 2*f(a+i*h);	//xk = a + kh ,
    }
    ans *= h/2;
    return ans;
}

复化辛普森公式

记字段 $x_k,x_{k+1}]$ 的中点为 $x_{k+\frac 12}=\frac {x_k+x_{k+1}}2 =a+\frac {(2k+1)h}2$ ,
在每个子区间 $x_k,x_{k+1}]$ ，用辛普森（Simpson）公式

$\int_{a}^b f(x)dx\approx \frac {b-a}6[f(a)+4f(\frac {a+b}2)+f(b)]$

得

$S_n = \sum_{k=0}^{n-1} \frac h6 [f(x_k)+4f(x_{k+ \frac 12})+f(x_{k+1})]$

按上文的方式展开化简得

$S_n=\frac h6[f(a)+4\sum_{k=0}^{n-1} f(x_{k+ \frac 12})+2\sum_{k=1}^{n-1}f(x_k)+f(b)]$

令
$\sum_{k=0}^{n-1} f(x_{k+ \frac 12})=\sum_{k=1}^{n} f(x_{k- \frac 12})\\ 2\sum_{k=1}^{n-1}f(x_k)=2\sum_{k=1}^{n}f(x_k)-2f(x_n)=2\sum_{k=1}^{n}f(x_k)-2f(b)\\$

最后整理得

$S_n=\frac h6 \{f(a)-f(b)+\sum_{k=1}^{n}[ 4f(x_{k- \frac 12})+2f(x_k)] \}$

易知 $x_{k-\frac 12} =a+\frac {(2k-1)h}2$ ，复化辛普森公式的实现与复化梯形公式类似

代码如下：

double CompositeSimpson(double a,double b,double n,double(*f)(double x))
{
    double h = (b - a) / n;    
    double ans = f(a) - f(b);
    double x = a;
    for(int i=1;i<=n;i++)
    {
        ans += 4*f(a+(2*i-1)/2.0*h);	// x(k-1/2)
        ans += 2*f(a+i*h);				// xk
    }
    ans *= h/6;
    return ans;
}

测试

用复化梯形公式、复化辛普森公式求积分 $\int_1^2 \frac 23 x^3 e^{x^2}dx=e^4$ ，
取n=10

完整代码

#include <iostream>
#include <cmath>
using namespace std;

// 复化梯形公式
double CompositeTrapezoidal(double a,double b,double n,double(*f)(double x))
{
    double h = (b - a) / n;    
    double ans = f(a) + f(b);
    for(int i = 1;i < n;i++)
    {
        ans += 2*f(a+i*h);
    }
    ans *= h/2;
    return ans;
} 
// 复化辛普森公式
double CompositeSimpson(double a,double b,double n,double(*f)(double x))
{
    double h = (b - a) / n;    
    double ans = f(a) - f(b);
    double x = a;
    for(int i=1;i<=n;i++)
    {
        ans += 4*f(a+(2*i-1)/2.0*h);
        ans += 2*f(a+i*h);
    }
    ans *= h/6;
    return ans;
}

int main()
{
    double a,b,h;

    int n;
     cout<<"输入a :";cin>>a;
     cout<<"输入b :";cin>>b;
     cout<<"输入n :";cin>>n;
    
    auto f=[](double x){return (2.0/3*x*x*x*exp(x*x));};
    cout<<"精确值"<< exp(4)<< endl;
    cout<<"使用复化梯形公式求得的结果:"<<CompositeTrapezoidal(a,b,n,f)
        <<",绝对误差为: "<<abs(CompositeTrapezoidal(a,b,n,f)-exp(4))<<endl;
    cout<<"使用复化辛普森公式求得的结果: "<<CompositeSimpson(a,b,n,f)
        <<",绝对误差为: "<<abs(CompositeSimpson(a,b,n,f)-exp(4))<<endl;
    return 0;
}