【老生谈算法】matlab实现自适应Simpson积分算法源码——积分算法

自适应Simpson积分算法（MATLAB及C++实现代码）

---

1、文档下载：

本算法已经整理成文档如下，有需要的朋友可以点击进行下载

---

2、算法详解：

对于类似问题，改一下积分函数和区间即可。

针对问题：数学上已经证明了

成立，所以可以通过数值积分来求的近似值。
试利用自适应Simpson算法计算积分近似值。

C++版：(直接复制粘贴在VC++6.0即可运行)

/*用自适应Simpson积分方法计算积分值*/
#include
#include
int n=0;       //设置全局变量n,用来记录最高迭代次数，避免递归一直进行下去。
double pi=3.141592653589793238462643 ;    //设置近似精确值，用以比较
double e1=0.00001  ; //设置误差容限为10^-5
double f(double);       //要积分的函数 
double Simpson (double,double,double,double);   // 迭代函数

using namespace std;

//主函数
 int main()
 {
	 double a=0,b=1,t,h,S;//积分区间
	
	 
	 h=(b-a)/2;
	 S=h/3*(f(a)+f(b)+4*f((a+b)/2));  //第一次Simpson公式积分值
	 
	 t=Simpson(a,b,e1,S);
	 cout<<"积分值为："<<t<<endl;
	 cout<<"最大迭代次数为："<<n<<endl;
	 cout<<"设置误差容限为"<<e1<<"\n误差为:"<<pi-t<<endl;

	 return 0; 
 }
 
 
 //子函数1（积分函数）
 double f(double x)
 {return 4/(1+x*x);}
 
 
 //子函数2（迭代函数）
 double Simpson (double A,double B,double e,double S)
 {double h,S1,S2;
   h=(B-A)/2;
   n++;    //统计迭代次数
   if(n>500)   {cout<<"方法有误，跳出递归"<<endl; return 0;}
 
 
   S1=h/6*(f(A)+f(A+h)+4*f(A+h/2));      // 在[A,(A+B)/2] 区间上计算Simpson积分值
   S2=h/6*(f(A+h)+f(B)+4*f(A+3/2*h));   // 在[(A+B)/2，B] 区间上计算Simpson积分值
   if(fabs(S-S1-S2)<15*e)                
	     return S1+S2;                  //如果满足误差容限要求 ，就以S1+S2作为此时对应区间上的函数的近似值
	else    
		return Simpson(A,(A+B)/2,e/2,S1)+Simpson((A+B)/2,B,e/2,S2); //递归调用
 }
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

MATLAB版： （两个函数文件加一个脚本文件）

1.编写积分函数文件：

function y =f(x)
y=4./(1+x.^2);
end
1
2
3

2.编写Simpson 迭代函数文件

function  y = Simpson( A,B,e,S )
h=(B-A)/2;
S1=h/6*(f(A)+f(A+h/2)+4*f(A+h/2));
S2=h/6*(f(A+h)+4*af(A+3/2*h)+f(B));
if  abs(S-S1-S2)<10*e    y= S1+S2;
else  y=Simpson(A,(A+B)/2,e/2,S1)+Simpson((A+B)/2,B,e/2,S2);
end
 
end
1
2
3
4
5
6
7
8
9

3.编写脚本调用文件

tic 
clear;    
a=0; b=1;     %积分区间
e=0.0000001;   %误差容限
h=(b-a)/2;       
S=h/3*(f(a)+f(b)+4*f(1/2*(a+b)));  %第一次Simpson积分值
t=Simpson(a,b,e,S)            %最终自适应方法积分值
abs(pi-t)                         %实际误差
e                             %设置的误差容限
toc   %返回所用时间
1
2
3
4
5
6
7
8
9
10

亲测可用。这两个代码本质上是一样的。我先用C++语言写好，然后又换用成MATLAB语言。

MATLAB好像可以把误差容限调到10^-7以下，而C++ 则只能到10^-5左右。原因不甚了解，猜测可能是由于C++计算时字节长度不够，导致精度不够，要递归调用很多次才能达到所需精度。