标签 算法 下的文章
算法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 1 条评论