标签 Simpson 下的文章
算法C/C++

算法(C++): 自适应辛普森（Simpson）积分和二重积分算法

概述

最近matlab中存在个二重积分的函数，携程c++的就给我整懵了，因为方法比较多，很容易结果对不上，之前的数值积分是用的龙贝格积分法，应该是前十位相同。这回辛普森的精确到9位，大概是没啥大问题。
辛普森积分法利用拉格朗日插值法来近似拟合原函数，然后通过对拟合函数的积分来简化近似。
近似公式如下：

$$ \int_{a}^{b} f(x) \mathrm{d} x \approx \frac{b-a}{6}\left(f(a)+4 f\left(\frac{a+b}{2}\right)+f(b)\right. $$

拉格朗日插值

拉格朗日插值利用二次曲线（抛物线）来拟合。
在求通过点集$(x_{1},y_{1})$,$(x_{2},y_{2})$,$(x_{3},y_{3})$,···,$(x_{n},y_{n})$的曲线上,拉格朗日利用多条二次曲线来拟合。首先二次曲线 $g = f(x)$ 满足

$$ f_{i}(x_{i}) \begin{cases} 1 \quad i = j\\ 0 \quad i \neq j\\ \end{cases} $$

因此$y_{n}f_{n}(x)$ 在$x_{n}$处，取值$y_{n}$，其余点取值为0。

拉格朗日通过对$f_{i}(x)$的构造如下：

$$ f_{i}(x)=\prod_{j=i}^{1 \leq j \leq n} \frac{\left(x-x_{j}\right)}{\left(x_{i}-x_{j}\right)} $$

最终得到了：

$$ f(x)=\sum_{i=1}^{n} y_{i}f_{i}(x) $$

辛普森公式

积分式：

$$ \int_{a}^{b} f(x)dx $$

第一我们需要用拉格朗日插值法对 $f(x)$ 进行近似计算，得出：

$$ f(x) \approx L_{n}(x) $$

$$ L_{n} = \sum_{i=1}^{n} y_{i}p_{i}(x) $$

$$ p_{i}(x)=\prod_{j=/i}^{1 \leq j \leq n} \frac{\left(x-x_{j}\right)}{\left(x_{i}-x_{j}\right)} $$

进而得出：

$$ f(x) \approx L_{n}(x) = \sum_{i=1}^{n} y_{i}p_{i}(x) $$

两边积分从而得出：

$$ \int_{a}^{b} f(x)dx \approx \int_{a}^{b} \sum_{i=1}^{n} y_{i}p_{i}(x) dx = \sum_{i=1}^{n} y_{i} \int_{a}^{b} p_{i}(x) dx $$

公式中的 $y_{i}=f(x_{i})& 。 令

$$ A_{i} = \int_{a}^{b} p_{i}(x) dx $$

那么就有：

$$ \int_{a}^{b} f(x)dx \approx \sum_{i=1}^{n} A_{i} f(x_{i}) $$

此外，如果在区间[a,b]中取n个点来拟合或者近似，则有 $h=\frac{b-a}{n}$,$X_{k}=a+kh$ , 所以：

$$ A_{i} = \int_{a}^{b} p_{i}(x) dx = h \int_{n}^{0} \prod_{j=1}^{n} \frac{(t-j)}{i-j}dt=(b-a)C_{a}^{(n)} $$

进而就有了：

$$ \int_{a}^{b} f(x)dx \approx \sum_{i=1}^{n} y_{i} A_{i} = (b-a)\sum_{i=0}^{n} C_{a}^{(n)} f(x_{i}) $$

该公式为牛顿-科特斯求积公式。
如果取3个点$x_{1}=a,x_{2}=\frac{a+b}{2},x_{3}=b$,则有下面：

$$ C_{0}^{(2)} = \frac{1}{6}, C_{1}^{(2)} = \frac{4}{6}, C_{2}^{(2)} = \frac{1}{6} $$

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

自适应辛普森积分

在辛普森公式中，三点的近似精度可能满足不了工程要求得精度，为此自适应辛普森积分会根据实际情况来自动调整精度。误差大的区域会划分多个区域。具体判断准则为：

$$ |S(a,c)+S(c,b)-S(a,b)| < 15 * eps $$

其中a , b为积分上下限，$c=\frac{a+b}{2}$，eps为精度

相关代码如下：

  • #include<iostream>

  • #include<stdlib.h>

  • using namespace std;

  • double Fun(double x){

  •     return x*x+x;//积分函数，自定义函数

  • }

  • double simpson(double a, double b){

  •    double c=(b+a)/2.0;

  •    return (Fun(a)+4*Fun(c)+Fun(b))*(b-a)/6.0;//求辛普森近似值

  • }

  • double adsp(double a,double b, double eps, double S){//自适应辛普森递归过程

  •     double c=(b+a)/2.0;

  •     double L=simpson(a,c),R=simpson(c,b);

  •     if(abs(L+R-S)<=15.0*eps){//判断准则

  •         return L+R+(L+R-S)/15.0;

  •     }

  •     return adsp(a,c,eps/2.0,L)+adsp(c,b,eps/2.0,R);

  • }

  • double intergation(double a, double b, double eps){//求积分

  •     return adsp(a,b,eps,simpson(a,b));

  • }

  • int main() {

  •     cout<<intergation(0, 1, 0.00001)<<endl;

  •     system("PAUSE");

  •     return 0;

  • }

二重积分

二重积分公式：

$$ \int_{a}^{b}\int_{c}^{d}f(x,y)dx dy = \int_{a}^{b}F_{x}(y)dy \approx \frac{b-1}{6}(F_{x}(a)+4F_{x}(\frac{a+b}{2})+F_{x}(b)) $$

$$ F_{x}(y) = \int_{d}^{c}f(x,y)dx \approx \frac{d-c}{6}(f(c,y)+4f(\frac{d+c}{2},y)+f(d,y)) $$

二重积分只要根据辛普森公式固定y的值后，对f ( x , y ) f(x,y)f(x,y)求关于x的积分,之后再对y求积分就ok了。

相关代码如下：
CommonAlgorithm.h

  • /**

  •  * @file CommonAlgorithm.h

  •  * @author ybw (root@bug-maker.com)

  •  * @brief 公共算法头文件

  •  * @version 0.1

  •  * @date 2022-01-20

  •  * 

  •  * @copyright Copyright (c) 2022 CIOMP

  •  * 

  •  */

  • #ifndef _COMMON_ALGORITHM_H_

  • #define _COMMON_ALGORITHM_H_

  • 

  • #include <iostream>

  • #include <string>

  • 

  • #ifndef M_PI

  • #define M_PI 3.1415926535897932384626433832795

  • #endif

  • 

  • typedef double(*cbp)(double);

  • typedef double(*cbp2)(double,double);

  • 

  • namespace CommonAlgorithm {

  •     /**

  •      * @brief 二重积分函数（对matlab）

  •      * 

  •      * @param cbp2 回调函数

  •      * @param xa x上限

  •      * @param xb x下限

  •      * @param ya y上线

  •      * @param yb y下限

  •      * @param eps 步长

  •      * @return double 

  •      */

  •     double integral2(double (* F)(double x, double y), double xa, double xb,double ya, double yb,double eps = 1.0e-8);

  •     double simpsonX2(double (* F)(double x, double y), double a, double b,double y);

  •     double adspX2(double (* F)(double x, double y), double a,double b,double y, double eps, double S);

  •     double inte2(double (* F)(double x, double y), double a, double b,double y, double eps);

  •     double simpsonY2(double (* F)(double x, double y), double xa, double xb,double ya, double yb,double eps);

  •     double adspY2(double (* F)(double x, double y), double xa, double xb,double ya, double yb, double eps,double S);

  •     double intergation2(double (* F)(double x, double y),double xa, double xb,double ya, double yb,double eps);

  • }

  • 

  • #endif

CommonAlgorithm.cpp

  • /**

  •  * @file CommonAlgorithm.cpp

  •  * @author ybw (root@bug-maker.com)

  •  * @brief 公共算法源文件

  •  * @version 0.1

  •  * @date 2022-01-20

  •  * 

  •  * @copyright Copyright (c) 2022 CIOMP

  •  * 

  •  */

  • #include "CommonAlgorithm.h"

  • 

  • double CommonAlgorithm::simpsonX2(double (* vFunction)(double x, double y), double a, double b,double y){//a,d为积分下、上限，y为被固定的y值

  •    double c=(b+a)/2.0;

  •    return (vFunction(a,y)+4*vFunction(c,y)+vFunction(b,y))*(b-a)/6.0;//对x的辛普森近似

  • }

  • 

  • double CommonAlgorithm::adspX2(double (* vFunction)(double x, double y), double a,double b,double y, double eps, double S){//对x的自适应辛普森递归

  •     double c = (b + a) / 2.0;

  •     double L = CommonAlgorithm::simpsonX2(vFunction,a,c,y),R=CommonAlgorithm::simpsonX2(vFunction,c,b,y);

  •     if(abs(L+R-S)<=15.0*eps){

  •         return L+R+(L+R-S)/15.0;

  •     }

  •     return CommonAlgorithm::adspX2(vFunction,a,c,y,eps/2.0,L)+CommonAlgorithm::adspX2(vFunction,c,b,y,eps/2.0,R);

  • }

  • double CommonAlgorithm::inte2(double (* vFunction)(double x, double y), double a, double b,double y, double eps){//固定y后，对x的积分

  •     return CommonAlgorithm::adspX2(vFunction,a,b,y,eps,CommonAlgorithm::simpsonX2(vFunction,a,b,y));

  • }

  • double CommonAlgorithm::simpsonY2(double (* vFunction)(double x, double y), double xa, double xb,double ya, double yb,double eps){

  •     double yc=(ya+yb)/2.0;

  •     return (CommonAlgorithm::inte2(vFunction,xa,xb,ya,eps)+4*CommonAlgorithm::inte2(vFunction,xa,xb,yc,eps)+CommonAlgorithm::inte2(vFunction,xa,xb,yb,eps))*(yb-ya)/6.0;//对y的辛普森近似

  • }

  • double CommonAlgorithm::adspY2(double (* vFunction)(double x, double y), double xa, double xb,double ya, double yb, double eps,double S){//对y的自适应辛普森递归

  •     double L=CommonAlgorithm::simpsonY2(vFunction,xa,xb,ya,(ya+yb)/2,eps);

  •     double R=CommonAlgorithm::simpsonY2(vFunction,xa,xb,(ya+yb)/2,yb,eps);

  •     if(abs(L+R-S)<=15.0*eps){

  •         return L+R+(L+R-S)/15.0;

  •     }

  •     return CommonAlgorithm::adspY2(vFunction,xa,xb,ya,(ya+yb)/2.0,eps/2.0,L)+CommonAlgorithm::adspY2(vFunction,xa,xb,(ya+yb)/2.0,yb,eps/2.0,R);

  • }

  • double CommonAlgorithm::intergation2(double (* vFunction)(double x, double y),double xa, double xb,double ya, double yb,double eps){//求二重积分

  •     return CommonAlgorithm::adspY2(vFunction,xa,xb,ya,yb,eps,CommonAlgorithm::simpsonY2(vFunction,xa,xb,ya,yb,eps));

  • }

  • 

  • double CommonAlgorithm::integral2(double (* vFunction)(double x, double y), double xa, double xb,double ya, double yb,double eps){

  •     return CommonAlgorithm::adspY2(vFunction,xa,xb,ya,yb,eps,CommonAlgorithm::simpsonY2(vFunction,xa,xb,ya,yb,eps));

  • }
阅读全文
2022-01-20 0 条评论