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