编程MATLAB程序实现复化梯形和辛普森数值积分
MATLAB数值分析实验二(复合梯形、辛普森和龙贝格求积,以及二重积分计算等)
佛山科学技术学院实验报告课程名称_______________ 数值分析________________________实验项目_______________ 数值积分____________________专业班级机械工程姓名余红杰学号2111505010 指导教师陈剑成绩日期月日一、实验目的b1、理解如何在计算机上使用数值方法计算定积分 a f ""X的近似值;2、学会复合梯形、复合Simpson和龙贝格求积分公式的编程与应用。
3、探索二重积分.11 f (x, y)dxdy在矩形区域D = {( x, y) | a _ x _ b, c _ y _ d}的数值D积分方法。
二、实验要求(1)按照题目要求完成实验内容;(2)写出相应的Matlab程序;(3)给出实验结果(可以用表格展示实验结果);(4)分析和讨论实验结果并提出可能的优化实验。
(5)写出实验报告。
三、实验步骤1、用不同数值方法计算积xln xdx =-- 0 9(1)取不同的步长h,分别用复合梯形及复合辛普森求积计算积分,给出误差中关于h的函数,并与积分精确值比较两公式的精度。
(2)用龙贝格求积计算完成问题(1 )。
2、给出一种求矩形区域上二重积分的复化求积方法,然后计算二重积分..e"y dxdy,其中积分区域D二{0乞x岂1,0岂y乞1}。
1.%lnt_t.m复化梯形:function F = Int_t(x1,x2,n)%复化梯形求积公式% x1,x2为积分起点和中点%分为n个区间,没选用步长可以防止区间数为非整数。
%样点矩阵及其函数值:x = lin space(x1,x2 ,n+1);y = f(x);m = len gth(x);%本题中用Matlab计算端点位置函数值为NaN,故化为零: y(1) = 0;y(m) = 0;%算岀区间长度,步长h:h = (x2 -x1)/n;a = [1 2*o nes(1,m-2) 1];%计算估计的积分值:F = h/2*sum(a.*y);%f.mfun cti on y = f(x)y = sqrt(x).*log(x);%run 11.mclc,clear;%分为10个区间,步长0.1的积分值:F = In t_t(0,1,10);F10 = F%分为100个区间F = In t_t(0,1,100);F100 = F%误差计算W10 = abs((-4/9)-F10);W100 = abs((-4/9)-F100);W = [W10 W100]%复化辛普森:%l nt_s.mfun cti on F = In t_s(x1,x2 ,n)%复化梯形求积公式% x1,x2区间,分为n个区间。
用复合梯形公式和复合辛普森公式求函数积分
《数值分析》实验报告(模板)
学号********班级信科121姓名张凯茜
【实验课题】用复合梯形公式和复合辛普森公式求函数积分
【实验目标】
1.掌握复合梯形公式与复合辛普森公式的基本思想。掌握常用的数值积分方法(特别是梯形法、Simpson方法、Cotes公式、Romberg算法以及Gauss求积公式)的原理。
【附程序】
复合梯形公式
functionT=comptra(a,b,tol)
h=b-a;
k=0;
T=((f(a)+f(b))*h)/2;
P=T+1;
whileabs(P-T)>tol
P=T;
m=0; h=h/2;
fori=1:2^k
m=m+f(a+(2*i-1)*h);
end
T=0.5*P+m*h; k=k+1;
2.学会用matlab编程实现用复合梯形公式与复合辛普森公式求积分。
3.熟悉matlab软件的使用,通过实验体会常用数值积分方法的逐步精致化过程。
【理论概述与算法描述】
1.根据梯形公式 ,将区间【a,b】划分为n等份,分点x(k)=a+kh,h=(b-a)/n,k=0,1,2,3,……,在每个区间【x(k),x(k+1)】(k=0,1,2……n-1)上采用梯形公式,得
end
复合辛普森公式
functionS=comsinp(a,b,tol)
h=b-a;
k=1;
S=((f(a)+f(b)+4*f((a+b)/2))*h)/6;
P=S;
whileabs(P-S)>tol
P=S;
matlab利用复合梯形公式计算积分
matlab利用复合梯形公式计算积分复合梯形公式是一种常用的数值积分方法,用于近似计算定积分。
在本文中,我们将使用MATLAB编程语言来实现复合梯形公式,并计算给定函数的积分。
首先,我们需要了解复合梯形公式的原理。
复合梯形公式是通过将积分区间划分为多个小区间,并在每个小区间上使用梯形面积来近似计算定积分。
具体而言,对于一个函数f(x),我们将积分区间[a, b]划分为n个小区间,每个小区间的宽度为h=(b-a)/n。
然后,我们可以使用以下公式来计算定积分的近似值:∫[a, b] f(x) dx ≈ h/2 * (f(a) + 2*f(x1) +2*f(x2) + ... + 2*f(xn-1) + f(b))其中,x1, x2, ..., xn-1是每个小区间的中点。
接下来,我们将使用MATLAB编程语言来实现复合梯形公式。
首先,我们需要定义要计算积分的函数f(x),以及积分区间[a, b]和划分的小区间数n。
```matlab\nfunction result =composite_trapezoidal(f, a, b, n)\n h = (b - a)/ n;\n x = a:h:b;\n result = h/2 * (f(a) +2*sum(f(x(2:end-1))) + f(b));\nend\n```在上述代码中,我们首先计算小区间的宽度h,并生成一个包含所有小区间的向量x。
然后,我们使用MATLAB的sum函数来计算除首尾之外的所有小区间上函数值的和,并将其乘以h/2得到积分的近似值。
接下来,我们可以定义要计算积分的函数f(x)。
例如,我们可以计算函数f(x) = x^2在积分区间[0, 1]上的积分。
```matlab\nfunction y = f(x)\n y =x.^2;\nend\n```最后,我们可以调用composite_trapezoidal函数来计算定积分的近似值。
复合梯形公式、复合辛普森公式matlab
复合梯形公式、复合⾟普森公式matlab 1. ⽤1阶⾄4阶Newton-Cotes公式计算积分程序:function I = NewtonCotes(f,a,b,type)%syms t;t=findsym(sym(f));I=0;switch typecase 1,I=((b-a)/2)*(subs(sym(f),t,a)+subs(sym(f),t,b));case 2,I=((b-a)/6)*(subs(sym(f),t,a)+4*subs(sym(f),t,(a+b)/2)+...subs(sym(f),t,b));case 3,I=((b-a)/8)*(subs(sym(f),t,a)+3*subs(sym(f),t,(2*a+b)/3)+...3*subs(sym(f),t,(a+2*b)/3)+subs(sym(f),t,b));case 4,I=((b-a)/90)*(7*subs(sym(f),t,a)+...32*subs(sym(f),t,(3*a+b)/4)+...12*subs(sym(f),t,(a+b)/2)+...32*subs(sym(f),t,(a+3*b)/4)+7*subs(sym(f),t,b));case 5,I=((b-a)/288)*(19*subs(sym(f),t,a)+...75*subs(sym(f),t,(4*a+b)/5)+...50*subs(sym(f),t,(3*a+2*b)/5)+...50*subs(sym(f),t,(2*a+3*b)/5)+...75*subs(sym(f),t,(a+4*b)/5)+19*subs(sym(f),t,b));case 6,I=((b-a)/840)*(41*subs(sym(f),t,a)+...216*subs(sym(f),t,(5*a+b)/6)+...27*subs(sym(f),t,(2*a+b)/3)+...272*subs(sym(f),t,(a+b)/2)+...27*subs(sym(f),t,(a+2*b)/3)+...216*subs(sym(f),t,(a+5*b)/6)+...41*subs(sym(f),t,b));case 7,I=((b-a)/17280)*(751*subs(sym(f),t,a)+...3577*subs(sym(f),t,(6*a+b)/7)+...1323*subs(sym(f),t,(5*a+2*b)/7)+...2989*subs(sym(f),t,(3*a+4*b)/7)+...1323*subs(sym(f),t,(2*a+5*b)/7)+...3577*subs(sym(f),t,(a+6*b)/7)+751*subs(sym(f),t,b));endsyms xf=exp(-x).*sin(x);a=0;b=2*pi;I = NewtonCotes(f,a,b,1)N=1:I =N=2:I =N=3:I =(pi*((3*3^(1/2)*exp(-(2*pi)/3))/2 - (3*3^(1/2)*exp(-(4*pi)/3))/2))/4N=4:I =(pi*(32*exp(-pi/2) - 32*exp(-(3*pi)/2)))/452. 已知,因此可以通过数值积分计算的近似值。
工程数值计算matlab实验报告——辛普森数值微分数值积分
《工程数值计算》上机实验报告(第二次)学生姓名 *** 班级 ******* 学号 ***********任课教师 *** 上机时间 2019 年 11月 5 日,报告完成 2019 年 11月 6 日1.实验目的:做功计算重物在力F(x)的拖动下,从x的起始点到最后点的移动过程中,力的大小及角度都是变化的,数据如表所示,试完成以下任务。
任务 2.1:试计算所做的功(数值积分实验)(1)分别采用梯形公式、辛普森公式和柯特斯公式计算,并进行比较(参考程序2.1a)。
(2)先对数据点进行多项式回归,再用int函数进行积分计算,并比较(参考程序2.1b)。
任务 2.2:试计算力在位移方向的变化率(数值微分实验)(1)基于给定的数据点,利用三点插值求各点处的导数(参考程序2.2)。
2.计算方法:针对实验任务,结合课堂内容,说明解决方法,如:采用何种理论,列出相关公式,说明计算步骤,写出程序框图等;任务2.1:数值积分计算原理:梯形公式辛普森公式柯特斯公式(n=4,即五个节点时)流程图:任务2.2:数值微分计算原理(三点插值):已知三个节点、、上的函数值,则流程图:3.程序设计:根据前面提到的计算方法编写程序;写出程序代码,并结合计算方法对程序中关键步骤进行必要的文字说明;代码:任务2.1:数值积分%数据输入clear;xi=[0 5 10 15 20 25 30 35 40];a=xi(1);b=xi(end);h=b-a;Fi=[0;9;13;14;10.5;12;5;8;9];Th=[0.5;1.4;0.75;0.9;1.3;1.48;1.5;1.1;0.8];%(1)%求x方向分力大小for i=1:length(xi)P(i)=Fi(i)*cos(Th(i));end%梯形公式,两个求积节点T=[0.5;0.5];%(系数)I1=h*(T(1)*P(1)+T(2)*P(9));%(梯形公式)%辛普森公式,三个求积节点S=[1/6;2/3;1/6];%(辛普森公式系数)I2=h*(S(1)*P(1)+S(2)*P(5)+S(3)*P(9));%(辛普森公式)%柯特斯公式,五个求积节点C=[7/90;16/45;2/15;16/45;7/90];%(柯特斯公式系数)I4=h*(C(1)*P(1)+C(2)*P(3)+C(3)*P(5)+C(4)*P(7)+C(5)*P(9));%(柯特斯公式)%显示结果I=[I1;I2;I4]%(2)多项式回归,再用int函数积分计算Ya=polyfit(xi,P,2);%(2次多项式回归)Yb=polyfit(xi,P,3);%(3次多项式回归)Yc=polyfit(xi,P,4);%(4次多项式回归)%构造回归多项式函数形式syms xLa=Ya(1)*x^2+Ya(2)*x+Ya(3);%(构造2次多项式函数)Lb=Yb(1)*x^3+Yb(2)*x^2+Yb(3)*x+Yb(4);%(构造3次多项式函数)Lc=Yc(1)*x^4+Yc(2)*x^3+Yc(3)*x^2+Yc(4)*x+Yc(5);%(构造4次多项式函数)%积分计算Ia=int(La,a,b);Ib=int(Lb,a,b);Ic=int(Lc,a,b);%显示计算结果II=[vpa(Ia);vpa(Ib);vpa(Ic)]%(将分数转化为小数输出)任务2.2:数值微分%数据输入clear;xi=[0 5 10 15 20 25 30 35 40];%(间隔为空格)a=xi(1);b=xi(end);Fi=[0;9;13;14;10.5;12;5;8;9];Th=[0.5;1.4;0.75;0.9;1.3;1.48;1.5;1.1;0.8];n=length(xi);%主体程序for i=1:nyi(i)=Fi(i)*cos(Th(i));endR=zeros(n,2);%(建立一个9*2的矩阵)%矩阵R的第一列记为点的序数for i=1:nR(i,1)=i;end%三点插值求导,导数值记入R矩阵的第二列中for i=1:nif i==1R(i,2)=(-3*yi(i)+4*yi(i+1)-yi(i+2))/(xi(i+2)-xi(i));%(第一个点的导数) else if i==nR(i,2)=(yi(i-2)-4*yi(i-1)+3*yi(i))/(xi(i)-xi(i-2));%(最后一个点的导数)elseR(i,2)=(yi(i+1)-yi(i-1))/(xi(i+1)-xi(i-1));%(中间点的导数)endendend%输出结果R4.结果分析:给出计算结果(可用数值、图表、曲线表示),并进行分析;任务2.1:数值积分由输出结果可以看出,梯形公式、辛普森公式、柯特斯公式的计算结果相差较大而2、3、4次多项式回归后再积分的结果相差较小,其中2、3次回归结果相差几乎可以不计。
matlab利用复合梯形公式计算积分
matlab利用复合梯形公式计算积分(原创版)目录1.MATLAB 简介2.复合梯形公式3.利用 MATLAB 计算积分的步骤4.举例说明5.总结正文1.MATLAB 简介MATLAB(Matrix Laboratory)是一款广泛应用于科学计算、数据分析、可视化等领域的软件。
它具有强大的数值计算和符号计算功能,使得用户可以方便地处理和解决各种数学问题。
在工程技术、自然科学等诸多领域,MATLAB 已成为研究和开发的重要工具。
2.复合梯形公式在 MATLAB 中,可以利用复合梯形公式(Composite Simpson Rule)来计算定积分。
复合梯形公式是一种高精度的数值积分方法,它将积分区间划分为多个子区间,然后在每个子区间上使用梯形公式进行近似计算。
最后将各子区间的结果加权求和,得到最终的积分结果。
3.利用 MATLAB 计算积分的步骤以下是使用 MATLAB 计算积分的具体步骤:(1)首先,需要导入 MATLAB 的符号计算工具箱,使用命令"syms"创建符号变量。
(2)定义被积函数。
例如,如果要计算积分∫(0, π) sin(x) dx,可以将 sin(x) 作为被积函数。
(3)使用"int"函数创建积分对象。
在 MATLAB 中,可以直接使用"int(sin(x), 0, pi)"创建积分对象。
(4)使用"quad"函数计算积分结果。
例如,可以输入"quad(int(sin(x), 0, pi))"来计算积分结果。
4.举例说明下面,我们通过一个具体的例子来说明如何使用 MATLAB 计算积分。
假设我们要计算定积分:∫(0, π) sin(x) dx。
首先,导入符号计算工具箱:"syms";然后,定义被积函数:"f(x) = sin(x)";接着,创建积分对象:"int(f(x), 0, pi)";最后,计算积分结果:"quad(int(f(x), 0, pi))"。
用matlab求数值积分的方法
用matlab求数值积分的方法
数值积分是一种求解定积分近似值的方法。
在实际应用中,很多复杂函数难以通过解析方法求得定积分,因此需要借助数值积分方法来求解。
Matlab作为一种高效的数值计算软件,提供了多种数值积分方法,包括梯形法、辛普森法、高斯积分法等。
下面分别介绍这些方法的具体实现。
梯形法:将积分区间等分成若干个小区间,每个小区间内的积分近似用其两端点的函数值的平均值。
最终将所有小区间的积分结果相加即为整个积分的近似值。
辛普森法:同样将积分区间等分成若干个小区间,每三个小区间内的积分近似用一个二次函数来拟合。
最终将所有小区间的积分结果相加即为整个积分的近似值。
高斯积分法:通过将积分区间变换为[-1,1]区间上的积分,利用预先计算好的高斯积分点和权重,将原函数在[-1,1]区间上积分近似为高斯点的函数值和权重之积的加权和。
以上就是Matlab中求解数值积分的三种常用方法。
不同方法在精度和计算速度上也有所差别,具体使用时可以根据实际需求进行选择。
- 1 -。
用matlab求数值积分的方法
用matlab求数值积分的方法
数值积分也称为数值积分法,是一种用计算机来近似求解定积分的方法。
在MATLAB中,可以使用三种数值积分方法:梯形法、辛普森法和积分变换法。
梯形法是最简单的数值积分方法之一,它通过将被积函数在区间上近似为一条直线,然后计算这条直线下的面积来近似定积分。
在MATLAB中,可以使用trapz函数来使用梯形法进行数值积分。
辛普森法是梯形法的改进版,它通过将被积函数在区间上近似为一个二次函数,然后计算这个二次函数下的面积来近似定积分。
在MATLAB中,可以使用quad函数来使用辛普森法进行数值积分。
积分变换法是一种更加精确的数值积分方法,它通过将被积函数进行一定的变换,然后将变换后的函数在区间上近似为一个多项式函数,最后计算这个多项式函数下的面积来近似定积分。
在MATLAB中,可以使用quadgk函数来使用积分变换法进行数值积分。
总之,MATLAB提供了许多数值积分的函数,选择合适的数值积分方法可以根据具体问题的要求来确定。
利用复化梯形公式、复化simpson公式计算积分
1、 利用复化梯形公式、复化simpson 公式计算积分2、 比较计算误差与实际误差取n=2,3,…,10分别利用复化梯形公式、复化simpson 公式计算积分1 I = x 2dx ,并与真值进行比较,并画出计算误差与实际误差之间的曲线。
利用复化梯形公式的程序代码如下: function f=fx(x) f=x.A 2;R=on es(1,9)*(-(b-a)/12*h.A 2*2); %积分余项(计算误差)true=quad(@fx,0,1); %积分的真实值A=T-true; %计算的值与真实值之差(实际误差)x=li nspace(0,1,9);plot(x,A,'r',x,R,'*')%将计算误差与实际误差用图像画出来 注:由于被积函数是x.A2,它的二阶倒数为 2,所以积分余项为:(-(b-a)/12*h.A 2*2) 实 验 原 理 ( a=0; b=1; T=[]; for n=2:10; %积分下线 %积分上线 %用来装不同n 值所计算出的结果 h=(b-a)/n; %步长 x=zeros(1, n+1); for i=1: n+1 x(i)=a+(i-1)*h; end y=x.A2; t=0; for i=1: n t=t+h/2*(y(i)+y(i+1)); end T=[T,t]; end%给节点定初值 %给节点赋值 %给相应节点处的函数值赋值 %利用复化梯形公式求值 %把不同n 值所计算出的结果装入 T 中 实验目的或 %首先建立被积函数,以便于计算真实值。
2法二:a=0;b=1;T=[];for n=2:10h=(b-a)/(2* n); x=zeros(1,2* n+1);for i=1:2* n+1x(i)=a+(i-1)*h;endy=x.A4;t=y(1)+y(2* n+1);for i=1: nt=t+4*y(2*i)+2*y(2*i-1);endT=[T,h/3*t];endtrue=quad(@fx1,0,1);A=T-true;x=li nspace(0,1,9);plot(x,A)此法与第一种一样,只是所用的表达式不同。
matlab编程积分复合辛普森公式
matlab编程积分复合辛普森公式编程求解复合辛普森公式是一种常用的数值积分方法,在MATLAB中可以通过编写相应的代码实现。
本文将详细介绍如何使用MATLAB编程计算复合辛普森公式,并给出具体的代码实现。
我们需要了解什么是复合辛普森公式。
复合辛普森公式是一种数值积分方法,用于近似计算函数的定积分。
它是在区间[a, b]上使用多个小区间进行逼近,而不是直接在整个区间上进行逼近。
这种方法的优势在于可以提高计算精度,并且对于复杂的函数也能够得到较好的近似结果。
我们需要将整个区间[a, b]划分为n个小区间。
每个小区间的长度为h=(b-a)/n。
然后,我们可以使用复合辛普森公式来近似计算每个小区间上的定积分。
复合辛普森公式的表达式为:I = (h/6)*(f(a) + 4*f((a+b)/2) + f(b))其中,f(x)是要求解的函数。
根据复合辛普森公式的定义,我们需要对每个小区间应用该公式进行求解,并将结果累加得到最终的积分值。
在MATLAB中,我们可以通过编写以下代码来实现复合辛普森公式的求解:```matlabfunction I = composite_simpson(f, a, b, n)h = (b-a)/n;x = a:h:b;y = f(x);I = 0;for i = 1:nI = I + h/6*(y(i) + 4*y(i+1) + y(i+2));endend```上述代码中,函数composite_simpson接受四个参数:函数f、积分区间的起点a、终点b和划分的小区间数n。
其中,函数f是一个函数句柄,表示要求解的函数。
在代码中,我们首先计算出每个小区间的长度h,并生成对应的x值。
然后,通过调用函数f计算出对应的y值。
接下来,我们使用循环对每个小区间应用复合辛普森公式,并将结果累加到变量I中。
最后,我们将得到的积分值I作为函数的输出。
在使用该函数时,我们需要先定义要求解的函数,并将其作为参数传递给composite_simpson函数。
复化梯形公式matlab
复化梯形公式matlab在MATLAB中,可以使用复化梯形公式来进行数值积分的计算。
复化梯形公式是将积分区间分割成多个小区间,在每个小区间上采用梯形面积逼近曲线下的积分值。
以下是使用MATLAB编写的复化梯形公式的示例代码:```matlabfunction result = composite_trapezoidal(f, a, b, n)h = (b - a) / n; % 计算每个小区间的宽度% 计算每个小区间的积分值,并将其累加得到最终结果result = 0;for i = 1:nx_i = a + (i-1) * h; % 当前小区间的起点x_j = a + i * h; % 当前小区间的终点% 使用梯形公式计算当前小区间的积分值integral_i = (f(x_i) + f(x_j)) * h / 2;% 将当前小区间的积分值累加到总结果中result = result + integral_i;endend```在上述代码中,`f` 是要计算积分的函数,`a` 和 `b` 是积分区间的起点和终点,`n` 是将积分区间划分成的小区间数目。
你可以根据实际需求,将自己的函数替换到 `f` 的位置,并调用 `composite_trapezoidal` 函数来计算数值积分的近似值。
例如,假设要计算函数 `f(x) = x^2` 在区间 `[0, 1]` 上的数值积分,可以使用以下代码进行计算:```matlabf = @(x) x^2; % 定义要计算积分的函数a = 0; % 积分区间起点b = 1; % 积分区间终点n = 100; % 将积分区间划分为100个小区间result = composite_trapezoidal(f, a, b, n); % 使用复化梯形公式计算积分近似值disp(result); % 显示计算结果```运行上述代码,就可以得到函数 `f(x) = x^2` 在区间 `[0, 1]` 上的数值积分的近似值。
MATLAB实现复化梯形公式复化SIMPSON公式以及ROMBERG积分
MATLAB实现复化梯形公式复化SIMPSON公式以及ROMBERG积分复化梯形公式、复化SIMPSON公式和ROMBERG积分是常用的数值积分方法,用于对定积分进行数值近似计算。
下面将介绍MATLAB实现这三种方法的具体步骤。
复化梯形公式使用多个等距的子区间进行近似计算,然后将子区间上的梯形面积求和。
MATLAB代码如下:```matlabh=(b-a)/n;%子区间宽度x=a:h:b;%子区间节点y=f(x);%子区间节点对应的函数值result = h * (sum(y) - (y(1) + y(end)) / 2); % 计算近似积分值end```复化SIMPSON公式同样使用多个等距的子区间进行近似计算,但是每个子区间上使用二次多项式拟合。
MATLAB代码如下:```matlabh=(b-a)/n;%子区间宽度x=a:h:b;%子区间节点y=f(x);%子区间节点对应的函数值result = (h / 3) * (y(1) + y(end) + 4 * sum(y(2:2:end-1)) + 2 * sum(y(3:2:end-2))); % 计算近似积分值end```3. ROMBERG积分(Romberg Integration)ROMBERG积分是一种逐次精化的数值积分方法,通过不断提高梯形法则的阶数进行近似计算。
MATLAB代码如下:```matlabfunction result = romberg_integration(f, a, b, n)R = zeros(n, n); % 创建一个n*n的矩阵用于存储结果h=b-a;%子区间宽度R(1,1)=(h/2)*(f(a)+f(b));%计算初始近似积分值for j = 2:nh=h/2;%缩小子区间宽度sum = 0;for i = 1:2^(j-2)sum = sum + f(a + (2 * i - 1) * h);endR(j, 1) = 0.5 * R(j-1, 1) + (h * sum); % 使用梯形法则计算积分值for k = 2:jR(j, k) = R(j, k-1) + (R(j, k-1) - R(j-1, k-1)) / ((4^k) - 1); % 使用Romberg公式计算积分值endendresult = R(n, n); % 返回最终近似积分值end```以上是MATLAB实现复化梯形公式、复化SIMPSON公式以及ROMBERG积分的代码。