MATLAB实验三-定积分的近似计算

matlab蒙特卡洛法求定积分摘要：一、蒙特卡洛法简介二、蒙特卡洛法求定积分的步骤1.确定积分区间2.生成随机点3.计算原函数在随机点的值4.计算积分三、Matlab 编程实现蒙特卡洛法求定积分1.编写蒙特卡洛积分函数2.生成随机点3.计算原函数在随机点的值4.计算积分四、结论正文：蒙特卡洛法是一种通过随机抽样来估算数学期望值的方法。

在求解定积分问题时，蒙特卡洛法可以通过模拟原函数在区间上的取值，用随机点代替实际点，从而近似计算定积分。

Matlab 提供了丰富的工具箱，可以方便地实现蒙特卡洛法求定积分。

首先，我们需要确定积分区间。

例如，对于积分区间[a, b]，我们需要确定a 和b 的值。

接下来，生成随机点。

在Matlab 中，我们可以使用rand 函数生成随机数。

在[a, b] 区间内生成足够多的随机点，例如n 个点。

然后，计算原函数在随机点的值。

这需要我们先编写一个原函数的Matlab 函数。

例如，如果原函数是f(x)，我们可以使用f(x) = ...来计算原函数在随机点的值。

最后，计算积分。

根据蒙特卡洛法的原理，我们可以用生成随机点的个数n 来近似计算定积分。

具体地，定积分的值约等于(1/n) * Σ[f(xi)]，其中xi 是随机点。

在Matlab 中，我们可以使用sum 函数计算Σ[f(xi)]。

综上所述，我们可以通过以下步骤在Matlab 中实现蒙特卡洛法求定积分：1.确定积分区间。

2.生成随机点。

3.计算原函数在随机点的值。

4.计算积分。

通过这种方法，我们可以有效地求解定积分问题。

实验三定积分的近似计算一、问题背景与实验目的利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数表达出来的情形．如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法．在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分．本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法．对于定积分的近似数值计算，Matlab有专门函数可用．二、相关函数（命令）及简介1．sum(a)：求数组a的和．2．format long：长格式，即屏幕显示15位有效数字．（注：由于本实验要比较近似解法和精确求解间的误差，需要更高的精度）．3．double()：若输入的是字符则转化为相应的ASCII码；若输入的是整型数值则转化为相应的实型数值．4．quad()：抛物线法求数值积分．格式：quad(fun,a,b) ，注意此处的fun是函数，并且为数值形式的，所以使用*、/、^等运算时要在其前加上小数点，即.*、./、.^等．例：Q = quad('1./(x.^3-2*x-5)',0,2);5．trapz()：梯形法求数值积分．格式：trapz(x,y)其中x为带有步长的积分区间；y为数值形式的运算（相当于上面介绍的函数fun）例：计算0sin()dx xπ⎰x=0:pi/100:pi;y=sin(x);trapz(x,y)6．dblquad()：抛物线法求二重数值积分．格式：dblquad(fun,xmin,xmax,ymin,ymax)，fun可以用inline定义，也可以通过某个函数文件的句柄传递．例1：Q1 = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)顺便计算下面的Q2，通过计算，比较Q1 与Q2结果（或加上手工验算），找出积分变量x、y的上下限的函数代入方法．Q2 = dblquad(inline('y*sin(x)'), 0, pi, pi, 2*pi)例2：Q3 = dblquad(@integrnd, pi, 2*pi, 0, pi)这时必须存在一个函数文件integrnd.m：function z = integrnd(x, y) z = y*sin(x);7．fprintf （文件地址，格式，写入的变量）：把数据写入指定文件．例：x = 0:.1:1; y = [x; exp(x)];fid = fopen('exp.txt','w'); %打开文件 fprintf(fid,'%6.2f %12.8f\n',y); %写入 fclose(fid) %关闭文件 8．syms 变量1 变量2 …：定义变量为符号． 9．sym('表达式')：将表达式定义为符号．解释：Matlab 中的符号运算事实上是借用了Maple 的软件包，所以当在Matlab 中要对符号进行运算时，必须先把要用到的变量定义为符号． 10．int(f,v,a,b)：求f 关于v 积分，积分区间由a 到b ．11．subs(f ，'x'，a)：将 a 的值赋给符号表达式 f 中的 x ，并计算出值．若简单地使用subs(f)，则将f 的所有符号变量用可能的数值代入，并计算出值．三、实验内容1． 矩形法根据定积分的定义，每一个积分和都可以看作是定积分的一个近似值，即1()d ()nbi i ai f x x f x ς==∆∑⎰在几何意义上，这是用一系列小矩形面积近似小曲边梯形的结果，所以把这个近似计算方法称为矩形法．不过，只有当积分区间被分割得很细时，矩形法才有一定的精确度．针对不同i ς的取法，计算结果会有不同，我们以 120d 1xx +⎰为例（取100=n ），（1） 左点法：对等分区间b x i n ab a x x a x n i =<<-+=<<<=ΛΛ10，在区间],[1i i x x -上取左端点，即取1-=i i x ς，12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78789399673078， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差0.7878939967307840.0031784ππ-=≈（2）右点法：同（1）中划分区间，在区间],[1i i x x -上取右端点，即取i i x =ς，12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78289399673078， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差 0.7828939967307840.0031884ππ-=≈（3）中点法：同（1）中划分区间，在区间1[,]i i x x -上取中点，即取12i ii x x ς-+=， 12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78540024673078， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差 60.7854002467307842.653104ππ--=≈⨯如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数，那么可以期望得到比矩形法效果好得多的近似计算公式．下面介绍的梯形法和抛物线法就是这一指导思想的产物．2． 梯形法等分区间b x i n a b a x x a x n i =<<-+=<<<=ΛΛ10，nab x -=∆ 相应函数值为n y y y ,,,10Λ（n i x f y i i ,,1,0),(Λ==）．曲线)(x f y =上相应的点为n P P P ,,,10Λ（n i y x P i i i ,,1,0),,(Λ==）将曲线的每一段弧i i P P 1-用过点1-i P ，i P 的弦i i P P 1-（线性函数）来代替，这使得每个],[1i i x x -上的曲边梯形成为真正的梯形，其面积为x y y ii ∆⨯+-21，n i ,,2,1Λ=． 于是各个小梯形面积之和就是曲边梯形面积的近似值，11 11()d ()22nnbi i i i ai i y y x f x x x y y --==+∆≈⨯∆=+∑∑⎰， 即11 ()d ()22bn n ay y b a f x x y y n --≈++++⎰L ， 称此式为梯形公式．仍用 12 0d 1x x +⎰的近似计算为例，取100=n ，10112 0d ()122n n y y x b a y y x n --≈++++=+⎰L 0.78539399673078， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差 60.7853939967307845.305104ππ--=≈⨯很显然，这个误差要比简单的矩形左点法和右点法的计算误差小得多．3． 抛物线法由梯形法求近似值，当)(x f y =为凹曲线时，它就偏小；当)(x f y =为凸曲线时，它就偏大．若每段改用与它凸性相接近的抛物线来近似时，就可减少上述缺点，这就是抛物线法．将积分区间],[b a 作n 2等分，分点依次为b x i n a b a x x a x n i =<<-+=<<<=2102ΛΛ，nab x 2-=∆， 对应函数值为n y y y 210,,,Λ（n i x f y i i 2,,1,0),(Λ==），曲线上相应点为n P P P 210,,,Λ（n i y x P i i i 2,,1,0),,(Λ==）．现把区间],[20x x 上的曲线段)(x f y =用通过三点),(000y x P ，),(111y x P ，),(222y x P 的抛物线)(12x p x x y =++=γβα来近似代替，然后求函数)(1x p 从0x 到2x 的定积分：21 ()d x x p x x =⎰22 ()d x x x x x αβγ++=⎰)()(2)(30220223032x x x x x x -+-+-γβα]4)(2)()()[(62022022202002γβαγβαγβα++++++++++-=x x x x x x x x x x 由于2201x x x +=，代入上式整理后得 21 ()d x x p x x ⎰)](4)()[(612122202002γβαγβαγβα++++++++-=x x x x x x x x )4(621002y y y x x ++-=)4(6210y y y nab ++-= 同样也有422 ()d x x p x x ⎰)4(6432y y y n ab ++-=……222 ()d n n x nx p x x -⎰)4(621222n n n y y y nab ++-=-- 将这n 个积分相加即得原来所要计算的定积分的近似值：22222212 11()d ()d (4)6ii nnbx i i i i ax i i b af x x p x x y y y n---==-≈=++∑∑⎰⎰， 即021******* ()d [4()2()]6bn n n ab af x x y y y y y y y y n---≈++++++++⎰L L 这就是抛物线法公式，也称为辛卜生（Simpson ）公式．仍用 12 0d 1x x +⎰的近似计算为例，取100=n ，102132124222 0d [4()2()]16n n n x b ay y y y y y y y x n ---≈+++++++++⎰L L=0.78539816339745，理论值 12 0d 14x x π=+⎰，此时计算的相对误差 160.7853981633974542.827104ππ--=≈⨯4. 直接应用Matlab 命令计算结果（1） 数值计算 120d .1xx +⎰ 方法1：int('1/(1+x^2)','x',0,1) (符号求积分)方法2：quad('1./(1+x.^2)',0,1) （抛物线法求数值积分）方法3：x=0:0.001:1; y=1./(1+x.^2);trapz(x,y) （梯形法求数值积分） （2）数值计算 212 01d d x x y y -+⎰⎰方法1：int(int('x+y^2','y',-1,1),'x',0,2) （符号求积分）方法2：dblquad(inline('x+y^2'),0,2,-1,1) （抛物线法二重数值积分）四、自己动手1． 实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算 120d 1xx +⎰，取258=n ，并比较三种方法的精确程度．2． 分别用梯形法与抛物线法，计算 2 1d xx⎰，取120=n ．并尝试直接使用函数trapz()、quad()进行计算求解，比较结果的差异．3． 试计算定积分 0sin d xx x+∞⎰．（注意：可以运用trapz()、quad()或附录程序求解吗？为什么？）4． 将 120d 1xx +⎰的近似计算结果与Matlab 中各命令的计算结果相比较，试猜测Matlab 中的数值积分命令最可能采用了哪一种近似计算方法？并找出其他例子支持你的观点．5． 通过整个实验内容及练习，你能否作出一些理论上的小结，即针对什么类型的函数（具有某种单调特性或凹凸特性），用某种近似计算方法所得结果更接近于实际值？6． 学习fulu2sum.m 的程序设计方法，尝试用函数 sum 改写附录1和附录3的程序，避免for 循环．五、附录附录1：矩形法（左点法、右点法、中点法）（fulu1.m ） format long n=100;a=0;b=1;inum1=0;inum2=0;inum3=0; syms x fx fx=1/(1+x^2); for i=1:nxj=a+(i-1)*(b-a)/n; %左点 xi=a+i*(b-a)/n; %右点 fxj=subs(fx,'x',xj); %左点值fxi=subs(fx,'x',xi); %右点值fxij=subs(fx,'x',(xi+xj)/2); %中点值inum1=inum1+fxj*(b-a)/n;inum2=inum2+fxi*(b-a)/n;inum3=inum3+fxij*(b-a)/n;endinum1inum2inum3integrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum1 and real-value is about: %d\n\n',...abs((inum1-integrate)/integrate))fprintf('The relative error between inum2 and real-value is about: %d\n\n',...abs((inum2-integrate)/integrate))fprintf('The relative error between inum3 and real-value is about: %d\n\n',...abs((inum3-integrate)/integrate))附录2：梯形法（fulu2.m）format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x^2);for i=1:nxj=a+(i-1)*(b-a)/n;xi=a+i*(b-a)/n;fxj=subs(fx,'x',xj);fxi=subs(fx,'x',xi);inum=inum+(fxj+fxi)*(b-a)/(2*n);endinumintegrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))附录2sum：梯形法（fulu2sum.m），利用求和函数，避免for 循环format longn=100;a=0;b=1;syms x fxfx=1/(1+x^2);i=1:n;xj=a+(i-1)*(b-a)/n; %所有左点的数组xi=a+i*(b-a)/n; %所有右点的数组fxj=subs(fx,'x',xj); %所有左点值fxi=subs(fx,'x',xi); %所有右点值f=(fxi+fxj)/2*(b-a)/n; %梯形面积inum=sum(f) %加和梯形面积求解integrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))附录3：抛物线法（fulu3.m）format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x^2);for i=1:nxj=a+(i-1)*(b-a)/n; %左点xi=a+i*(b-a)/n; %右点xk=(xi+xj)/2; %中点fxj=subs(fx,'x',xj);fxi=subs(fx,'x',xi);fxk=subs(fx,'x',xk);inum=inum+(fxj+4*fxk+fxi)*(b-a)/(6*n);endinumintegrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))。

文章标题：探索matlab中的蒙特卡洛法求定积分在数学和计算科学中，求解定积分是一个常见的问题。

传统的数值积分方法中，蒙特卡洛法是一种非常有趣和强大的方法，能够对一些特殊的不易求解的定积分问题提供解决方案。

而在matlab这一强大的数学计算软件中，蒙特卡洛法同样有着广泛的应用。

1. 什么是蒙特卡洛法？蒙特卡洛法是一种基于随机采样的数值积分方法，其核心思想是利用随机抽样的方法逼近定积分的值。

具体来说，对于给定的函数$f(x)$以及区间$[a, b]$，蒙特卡洛法通过对函数在该区间上进行随机采样，并利用采样点的平均值来逼近定积分的值。

2. 在matlab中应用蒙特卡洛法在matlab中，可以利用蒙特卡洛法求解定积分问题。

通过生成服从均匀分布的随机数，并代入原函数，然后求解采样点的平均值，可以得到定积分的近似值。

matlab内置了丰富的数学计算和随机数生成函数，能够方便地实现蒙特卡洛法的计算。

3. 实例分析：使用matlab进行蒙特卡洛法求解定积分假设我们要求解函数$f(x)=x^2$在区间$[0, 1]$上的定积分，即$$\int_{0}^{1} x^2 \, dx$$我们可以在matlab中编写如下代码：```matlabN = 1000000; % 设定采样点的个数X = rand(1, N); % 生成均匀分布的随机数Y = X.^2; % 代入原函数integral_value = mean(Y); % 求解采样点的平均值```通过上述代码，我们得到了定积分的近似值integral_value。

在这个例子中，我们利用蒙特卡洛法求得了定积分的近似值。

4. 总结与展望通过本文的介绍，我们对matlab中蒙特卡洛法求解定积分的方法有了初步的了解。

蒙特卡洛法作为一种基于随机采样的数值积分方法，在matlab中有着广泛的应用。

在实际应用中，我们可以根据定积分的具体问题来灵活选择采样点的个数，并结合matlab强大的数学计算能力，在求解定积分问题中取得更加准确的结果。

用蒙特卡罗法估算定积分• 求解示例：程序是：clcn=10^6;maxfx=1;mimfx=0;a=0; b=1; m=0;d=1;c=[];for i=1:nx=a+rand*(b-a);y=rand*(maxfx-mimfx)+mimfx;if y<=x^2m=m+1;endendfprintf('积分值是：%f\n',m/n*(maxfx-mimfx)*(b-a))结果是： dx x 102积分值是：0.333116 积分值是：0.334410 积分值是：0.333132 积分值是：0.333022 积分值是：0.332899 积分值是：0.333604 积分值是：0.333535 积分值是：0.333117 积分值是：0.333312 积分值是：0.334002 积分值是：0.332966 积分值是：0.332274 积分值是：0.334320 积分值是：0.333143 积分值是：0.333144 积分值是：0.333840 积分值是：0.332857 积分值是：0.333487 积分值是：0.333332 积分值是：0.333501 积分值是：0.333358 积分值是：0.333129积分值是：0.333495 积分值是：0.334282 积分值是：0.333183 积分值是：0.333594 积分值是：0.333894 积分值是：0.333766 积分值是：0.333126 积分值是：0.332993 积分值是：0.333032 积分值是：0.333090 积分值是：0.332987 积分值是：0.333853 积分值是：0.333412 积分值是：0.332419 积分值是：0.333338 积分值是：0.333622 积分值是：0.334497 积分值是：0.333565 积分值是：0.333759 积分值是：0.332973 积分值是：0.332783 积分值是：0.332990积分值是：0.332452 积分值是：0.334463 积分值是：0.332927 积分值是：0.332749 积分值是：0.332269 积分值是：0.333389 积分值是：0.332691 积分值是：0.333038 积分值是：0.333918 积分值是：0.333517 积分值是：0.333182 积分值是：0.333919 积分值是：0.333463 积分值是：0.332823 积分值是：0.333679 积分值是：0.332826 >>Welcome To Download !!!欢迎您的下载，资料仅供参考！。

实验三定积分的近似计算一、问题背景与实验目的利用牛顿—莱布尼兹公式虽然可以精确地计算定积分的值，但它仅适用于被积函数的原函数能用初等函数表达出来的情形．如果这点办不到或者不容易办到，这就有必要考虑近似计算的方法．在定积分的很多应用问题中，被积函数甚至没有解析表达式，可能只是一条实验记录曲线，或者是一组离散的采样值，这时只能应用近似方法去计算相应的定积分．本实验将主要研究定积分的三种近似计算算法：矩形法、梯形法、抛物线法．对于定积分的近似数值计算，Matlab有专门函数可用．二、相关函数（命令）及简介1．sum(a)：求数组a的和．2．format long：长格式，即屏幕显示15位有效数字．（注：由于本实验要比较近似解法和精确求解间的误差，需要更高的精度）．3．double()：若输入的是字符则转化为相应的ASCII码；若输入的是整型数值则转化为相应的实型数值．4．quad()：抛物线法求数值积分．格式：quad(fun,a,b) ，注意此处的fun是函数，并且为数值形式的，所以使用*、/、^等运算时要在其前加上小数点，即.*、./、.^等．例：Q = quad('1./(x.^3-2*x-5)',0,2);5．trapz()：梯形法求数值积分．格式：trapz(x,y)其中x为带有步长的积分区间；y为数值形式的运算（相当于上面介绍的函数fun）例：计算0sin()dx xπ⎰x=0:pi/100:pi;y=sin(x);trapz(x,y)6．dblquad()：抛物线法求二重数值积分．格式：dblquad(fun,xmin,xmax,ymin,ymax)，fun可以用inline定义，也可以通过某个函数文件的句柄传递．例1：Q1 = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)顺便计算下面的Q2，通过计算，比较Q1 与Q2结果（或加上手工验算），找出积分变量x、y的上下限的函数代入方法．Q2 = dblquad(inline('y*sin(x)'), 0, pi, pi, 2*pi)例2：Q3 = dblquad(@integrnd, pi, 2*pi, 0, pi)这时必须存在一个函数文件integrnd.m：function z = integrnd(x, y) z = y*sin(x);7．fprintf （文件地址，格式，写入的变量）：把数据写入指定文件．例：x = 0:.1:1; y = [x; exp(x)];fid = fopen('exp.txt','w'); %打开文件 fprintf(fid,'%6.2f %12.8f\n',y); %写入 fclose(fid) %关闭文件 8．syms 变量1 变量2 …：定义变量为符号． 9．sym('表达式')：将表达式定义为符号．解释：Matlab 中的符号运算事实上是借用了Maple 的软件包，所以当在Matlab 中要对符号进行运算时，必须先把要用到的变量定义为符号． 10．int(f,v,a,b)：求f 关于v 积分，积分区间由a 到b ．11．subs(f ，'x'，a)：将 a 的值赋给符号表达式 f 中的 x ，并计算出值．若简单地使用subs(f)，则将f 的所有符号变量用可能的数值代入，并计算出值．三、实验内容1． 矩形法根据定积分的定义，每一个积分和都可以看作是定积分的一个近似值，即1()d ()nbi i ai f x x f x ς==∆∑⎰在几何意义上，这是用一系列小矩形面积近似小曲边梯形的结果，所以把这个近似计算方法称为矩形法．不过，只有当积分区间被分割得很细时，矩形法才有一定的精确度．针对不同i ς的取法，计算结果会有不同，我们以 120d 1xx +⎰为例（取100=n ），（1） 左点法：对等分区间b x i n ab a x x a x n i =<<-+=<<<=ΛΛ10，在区间],[1i i x x -上取左端点，即取1-=i i x ς，12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78789399673078， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差0.7878939967307840.0031784ππ-=≈（2）右点法：同（1）中划分区间，在区间],[1i i x x -上取右端点，即取i i x =ς，12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78289399673078， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差 0.7828939967307840.0031884ππ-=≈（3）中点法：同（1）中划分区间，在区间1[,]i i x x -上取中点，即取12i ii x x ς-+=， 12 01d ()1ni i i xf x x ς==∆≈+∑⎰0.78540024673078， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差 60.7854002467307842.653104ππ--=≈⨯如果在分割的每个小区间上采用一次或二次多项式来近似代替被积函数，那么可以期望得到比矩形法效果好得多的近似计算公式．下面介绍的梯形法和抛物线法就是这一指导思想的产物．2． 梯形法等分区间b x i n a b a x x a x n i =<<-+=<<<=ΛΛ10，nab x -=∆ 相应函数值为n y y y ,,,10Λ（n i x f y i i ,,1,0),(Λ==）．曲线)(x f y =上相应的点为n P P P ,,,10Λ（n i y x P i i i ,,1,0),,(Λ==）将曲线的每一段弧i i P P 1-用过点1-i P ，i P 的弦i i P P 1-（线性函数）来代替，这使得每个],[1i i x x -上的曲边梯形成为真正的梯形，其面积为x y y ii ∆⨯+-21，n i ,,2,1Λ=． 于是各个小梯形面积之和就是曲边梯形面积的近似值，11 11()d ()22nnbi i i i ai i y y x f x x x y y --==+∆≈⨯∆=+∑∑⎰， 即11 ()d ()22bn n ay y b a f x x y y n --≈++++⎰L ， 称此式为梯形公式．仍用 12 0d 1x x +⎰的近似计算为例，取100=n ，10112 0d ()122n n y y x b a y y x n --≈++++=+⎰L 0.78539399673078， 理论值 12 0d 14x x π=+⎰，此时计算的相对误差 60.7853939967307845.305104ππ--=≈⨯很显然，这个误差要比简单的矩形左点法和右点法的计算误差小得多．3． 抛物线法由梯形法求近似值，当)(x f y =为凹曲线时，它就偏小；当)(x f y =为凸曲线时，它就偏大．若每段改用与它凸性相接近的抛物线来近似时，就可减少上述缺点，这就是抛物线法．将积分区间],[b a 作n 2等分，分点依次为b x i n a b a x x a x n i =<<-+=<<<=2102ΛΛ，nab x 2-=∆， 对应函数值为n y y y 210,,,Λ（n i x f y i i 2,,1,0),(Λ==），曲线上相应点为n P P P 210,,,Λ（n i y x P i i i 2,,1,0),,(Λ==）．现把区间],[20x x 上的曲线段)(x f y =用通过三点),(000y x P ，),(111y x P ，),(222y x P 的抛物线)(12x p x x y =++=γβα来近似代替，然后求函数)(1x p 从0x 到2x 的定积分：21 ()d x x p x x =⎰22 ()d x x x x x αβγ++=⎰)()(2)(30220223032x x x x x x -+-+-γβα]4)(2)()()[(62022022202002γβαγβαγβα++++++++++-=x x x x x x x x x x 由于2201x x x +=，代入上式整理后得 21 ()d x x p x x ⎰)](4)()[(612122202002γβαγβαγβα++++++++-=x x x x x x x x )4(621002y y y x x ++-=)4(6210y y y nab ++-= 同样也有422 ()d x x p x x ⎰)4(6432y y y n ab ++-=……222 ()d n n x nx p x x -⎰)4(621222n n n y y y nab ++-=-- 将这n 个积分相加即得原来所要计算的定积分的近似值：22222212 11()d ()d (4)6ii nnbx i i i i ax i i b af x x p x x y y y n---==-≈=++∑∑⎰⎰， 即021******* ()d [4()2()]6bn n n ab af x x y y y y y y y y n---≈++++++++⎰L L 这就是抛物线法公式，也称为辛卜生（Simpson ）公式．仍用 12 0d 1x x +⎰的近似计算为例，取100=n ，102132124222 0d [4()2()]16n n n x b ay y y y y y y y x n ---≈+++++++++⎰L L=0.78539816339745，理论值 12 0d 14x x π=+⎰，此时计算的相对误差 160.7853981633974542.827104ππ--=≈⨯4. 直接应用Matlab 命令计算结果（1） 数值计算 120d .1xx +⎰ 方法1：int('1/(1+x^2)','x',0,1) (符号求积分)方法2：quad('1./(1+x.^2)',0,1) （抛物线法求数值积分）方法3：x=0:0.001:1; y=1./(1+x.^2);trapz(x,y) （梯形法求数值积分） （2）数值计算 212 01d d x x y y -+⎰⎰方法1：int(int('x+y^2','y',-1,1),'x',0,2) （符号求积分）方法2：dblquad(inline('x+y^2'),0,2,-1,1) （抛物线法二重数值积分）四、自己动手1． 实现实验内容中的例子，即分别采用矩形法、梯形法、抛物线法计算 120d 1xx +⎰，取258=n ，并比较三种方法的精确程度．2． 分别用梯形法与抛物线法，计算 2 1d xx⎰，取120=n ．并尝试直接使用函数trapz()、quad()进行计算求解，比较结果的差异．3． 试计算定积分 0sin d xx x+∞⎰．（注意：可以运用trapz()、quad()或附录程序求解吗？为什么？）4． 将 120d 1xx +⎰的近似计算结果与Matlab 中各命令的计算结果相比较，试猜测Matlab 中的数值积分命令最可能采用了哪一种近似计算方法？并找出其他例子支持你的观点．5． 通过整个实验内容及练习，你能否作出一些理论上的小结，即针对什么类型的函数（具有某种单调特性或凹凸特性），用某种近似计算方法所得结果更接近于实际值？6． 学习fulu2sum.m 的程序设计方法，尝试用函数 sum 改写附录1和附录3的程序，避免for 循环．五、附录附录1：矩形法（左点法、右点法、中点法）（fulu1.m ） format long n=100;a=0;b=1;inum1=0;inum2=0;inum3=0; syms x fx fx=1/(1+x^2); for i=1:nxj=a+(i-1)*(b-a)/n; %左点 xi=a+i*(b-a)/n; %右点 fxj=subs(fx,'x',xj); %左点值fxi=subs(fx,'x',xi); %右点值fxij=subs(fx,'x',(xi+xj)/2); %中点值inum1=inum1+fxj*(b-a)/n;inum2=inum2+fxi*(b-a)/n;inum3=inum3+fxij*(b-a)/n;endinum1inum2inum3integrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum1 and real-value is about: %d\n\n',...abs((inum1-integrate)/integrate))fprintf('The relative error between inum2 and real-value is about: %d\n\n',...abs((inum2-integrate)/integrate))fprintf('The relative error between inum3 and real-value is about: %d\n\n',...abs((inum3-integrate)/integrate))附录2：梯形法（fulu2.m）format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x^2);for i=1:nxj=a+(i-1)*(b-a)/n;xi=a+i*(b-a)/n;fxj=subs(fx,'x',xj);fxi=subs(fx,'x',xi);inum=inum+(fxj+fxi)*(b-a)/(2*n);endinumintegrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))附录2sum：梯形法（fulu2sum.m），利用求和函数，避免for 循环format longn=100;a=0;b=1;syms x fxfx=1/(1+x^2);i=1:n;xj=a+(i-1)*(b-a)/n; %所有左点的数组xi=a+i*(b-a)/n; %所有右点的数组fxj=subs(fx,'x',xj); %所有左点值fxi=subs(fx,'x',xi); %所有右点值f=(fxi+fxj)/2*(b-a)/n; %梯形面积inum=sum(f) %加和梯形面积求解integrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))附录3：抛物线法（fulu3.m）format longn=100;a=0;b=1;inum=0;syms x fxfx=1/(1+x^2);for i=1:nxj=a+(i-1)*(b-a)/n; %左点xi=a+i*(b-a)/n; %右点xk=(xi+xj)/2; %中点fxj=subs(fx,'x',xj);fxi=subs(fx,'x',xi);fxk=subs(fx,'x',xk);inum=inum+(fxj+4*fxk+fxi)*(b-a)/(6*n);endinumintegrate=int(fx,0,1)integrate=double(integrate)fprintf('The relative error between inum and real-value is about: %d\n\n',...abs((inum-integrate)/integrate))。