概述
最近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));
- }