《matlab教程加油.ppt》由会员分享，可在线阅读，更多相关《matlab教程加油.ppt（63页珍藏版）》请在三一办公上搜索。
1、实验4 积分计算,实验目的：,1.通过实验加深理解积分理论中分割、近似、求和、取极限的思想方法；,2.学习并掌握matlab求不定积分、定积分、二重积分、曲线积分的方法；,3.学习matlab命令sum，symsum，int。,学习 Matlab 命令计算不定积分定积分的概念不定积分和广义积分二重积分的计算曲线积分,实验内容：,1.学习Matlab命令,1).求和命令 sum 调用格式：,sum(x),给出向量x的各元素的累加和。若x为矩阵，则是一个元素为每列列和的行向量。,例1,x=1,2,3,4,5,6,7,8,9,10;sum(x),ans=55,例2,x=1,2,3;4,5,6;7,8
2、,9,x=1 2 3 4 5 6 7 8 9,sum(x),ans=12 15 18,symsum(s,n),symsum(s,k,m,n),当x的元素很有规律，例如为s(k)时，可用symsum求得x的各项和，n可取无穷。,2).求和命令 symsum 调用格式：,symsum(s(k),1,n)=s(1)+s(2)+.+s(n),symsum(s(k),k,m,n)=s(m)+s(m+1)+.+s(n),例3.,syms k nsymsum(k,1,10)symsum(k,2,10),ans=55ans=54,symsum(k2,k,1,n),ans=1/3*(n+1)3-1/2*(n+1
3、)2+1/6*n+1/6,例4.求下列部分和,syms n xs1=symsum(-1)(n+1)*x/(n*(n+2),n,1,30),s1=495/1984*x,syms n a ks2=symsum(-1)k*a*sin(k),k,0,n-1),syms ns3=symsum(1/n2,n,1,10),例5.讨论下列级数的敛散性,syms ns1=symsum(1/n2,n,1,inf),syms ns2=symsum(1/n,n,1,inf),syms n as3=symsum(an/n,n,1,inf),int(f(x),int(f(x,y),x),3).积分命令 int 调用格式：
4、,int(f(x),a,b),int(f(x,y),x,a,b),2.计算不定积分,例6,解：输入命令：,syms x;y=x2*log(x);y1=diff(y)%对函数求导y0=int(y1)%对导数积分的原函数y2=int(y),y2=1/3*x3*log(x)-1/9*x3,例7,解：输入命令：,syms x;syms a real;y=sqrt(a2-x2),(x-1)/(3*x-1)(1/3),x2*asin(x);int(y,x);,ans=1/2*x*(a2-x2)(1/2)+1/2*a2*asin(1/a2)(1/2)*x),1/15*(3*x-1)(5/3)-1/3*(3*
5、x-1)(2/3),1/3*x3*asin(x)+1/9*x2*(1-x2)(1/2)+2/9*(1-x2)(1/2),syms a x;f2=int(x3*(cos(a*x)2)f=simple(f2),1/16*(4*a3*x3*sin(2*a*x)+2*x4*a4+6*x2*a2*cos(2*a*x)-6*a*x*sin(2*a*x)+3-3*cos(2*a*x)/a4,f1=x4/8+(x3/(4*a)-3*x/(8*a3)*sin(2*a*x)+(3*x2/(8*a2)-3/(16*a4)*cos(2*a*x);r=simple(f-f1),r=3/16/a4,3 定积分的概念,x=
6、linspace(0,1,21);,定积分为一个和式极限，取f(x)=exp(x)，积分区间为0,1，等距划分为20个子区间,y=exp(x);,选取每个子区间的的端点，计算端点处的函数值,y1=y(1:20);s1=sum(y1)/20;,取区间的左端点乘以区间长度全部加起来,s1=1.6757,y2=y(2:21);s2=sum(y2)/20;,取区间的右端点乘以区间长度全部加起来,s2=1.7616,plot(x,y);hold onfor i=1:20fill(x(i),x(i+1),x(i+1),x(i),x(i),0,0,y(i),y(i),0,b)end,for i=1:20fi
7、ll(x(i),x(i+1),x(i+1),x(i),x(i),0,0,y(i+1),y(i+1),0,r)end,若取右端点，则,从图上可以看出：,当点取得越来越多时，s2-s1的值会越来越小，可试取50个点计算，看结果如何。下面按等分区间计算,syms k n;s=symsum(exp(k/n)/n,k,1,n);limit(s,n,inf),ans=exp(1)-1,结果与上面一样。,例9,解：,输入指令,syms x;I=int(exp(x),0,1),得结果:I=exp(1)-1,4 计算定积分和广义积分,例10,解：,输入指令,syms x;int(abs(x-1),0,2),an
8、s=1.,int()还可以求函数边界的定积分问题。,syms x t;f=(-2*(x2)+1)/(2*(x2)-3*x+1)2;I=simple(int(f,x,cos(t),exp(-2*t),I=-(exp(2*t)*cos(t)-1)*(-2*cos(t)+exp(2*t)/(-2+exp(2*t)/(-1+exp(2*t)/(2*cos(t)-1)/(cos(t)-1),例12,解：,对第一个积分输入指令：,syms x;syms p real;int(1/xp,x,1,inf),ans=limit(-(x-exp(p*log(x)/(p-1)/exp(p*log(x),x=inf)
9、,有结果看出，当p1时，ans=1/(p-1).,syms x;int(1/(x-1)2,0,2),ans=inf,对第二个积分，输入命令：,syms x;int(1/(2*pi)(1/2)*exp(-x2/2),-inf,inf),对第三个积分，输入命令：,ans=,例13,解：,输入指令：,syms x t;int(sin(x)/x,0,t),ans=sinint(t),help sinint SININT Sine integral function.SININT(x)=int(sin(t)/t,t,0,x).See also COSINT.Overloaded methods help
10、 sym/sinint.m,这类积分无法用初等函数或其值来表示。,5 数值积分,(1)梯形积分法-先将积分区间划分为几个小区间，用每个小区间上梯形面积之和作为积分的近似值。具体方法如下：,分点为：a=x0,x1,xn=b,将区间分成n个等长的小区间，区间长度为x=b-a/n，设函数y=f(x)对应的各分点的函数值为y0,y1,yn,每个小梯形的面积为：,从而有：,在编程时可采用如下方法估计误差：,可以证明梯形法最大误差为：,其中|f(x)|M,编程时可如下进行设误差为：,逐步计算A(n),若|A(n+1)-A(n)|,则以A(n)作为积分的近似值。,可先编写f.m文件定义函数y=f(x),然后
11、再编写A.m文件计算积分。,function A=A(t,e)A=sin(x)/x;for k=1:20 m=t;n=2k;x=linspace(0,t,n+1);y=sin(x)./x;h=1/n;T=h*(sum(y)-y(1)/2-y(n+1)/2);if(abs(m-t)=e)break;endend,z=trapz(x,y)x表示积分区间的离散化向量，y是与x同维数的被积函数向量，z返回积分的近似值.,例14.求积分,解：我们先用具有较大步长的积分区间的离散化向量求解,(高精度近似值为1.49364826),梯形积分法,clear;x=-1:0.5:1;y=exp(-x.2);z=t
12、rapz(x,y),（高精度近似值为1.49364826）,z=1.46274050,xx=-1:0.05:1;yy=exp(-x.2);z=trapz(xx,yy),z=1.49334167,plot(x,y,r-o,xx,yy),练习：用梯形法计算积分,（2）Simpson公式（抛物线法）-f 在a，b区间上积分的Simpson公式。,用通过三点(a,f(a),(a+b)/2,f(a+b)/2),(b,f(b)的抛物线围成的曲边梯形的面积来代替由 f 围成的曲边梯形的面积，由此获得积分的近似值.,-将区间a，b 分为n等分，步长为h=(ba)/n，分点为xk=a+kh(k=0,.,n)。在
13、每个小区间xk,xk+1上采用辛普森公式，得到小区间上 f 积分的近似值。对每个小区间上积分的近似值求和，得到f 在a，b区间积分的近似值，这就是复合辛普森公式：,可以证明最大误差不超过,其中M是|f(4)(x)|在a,b上的最大值。,(3)adaptive Simpson公式（自适应辛普生法）-通常被积函数f 在整个区间a，b 上变化不是很均匀的，如在某点附近函数变化非常急剧，而在其余地方上的变化比较平缓。这种情况用等距剖分小区间的复合求积公式不很合适。为了使积分达到预定的精度又要节省工作量，则可以在函数变化急剧的部分增多节点，即子区间分得细，而在函数变化平缓的地方减少节点，即子区间分得大，
14、这个方法称为自适应积分法(adaptive quadrature methods),-对于自适应辛普森方法，采用逐次将区间二等分的方法。为写法统一，将区间a，b记为a，a+h，其中h=ba为区间的长度，称为0级区间。在区间a,a+h上采用辛普森公式，把结果记作,将区间分成两个相等的子区间a,a+h/2和a+h/2,a+h,这两个子区间称为1级子区间，其长度为h/2.在每个子区间上 采用辛普森公式计算积分，然后相加并令,再将1级子区间中的一个或所有的两个二等分，所得的子区间称为2级子区间，其长度为h/22,.。如此继续下去，最后将区间 a,a+h分成n个子区间ai,ai+1(i=0,1,.,n
15、1)，这样有,子区间的长度一般是不同的，如果子区间 ai,ai+1是r级，则其长度为,实际上，区间的划分是根据被积函数的变化而定的，函数变化平缓的地方，子区间大，函数变化急剧的地方，子区间小。,设Sa,a+h表示积分I(f)的近似值，那么有,-自适应辛普森方法的matlab的命令为 z=quad(f,a,b,tol,trace),f(x)为被积函数a 为积分下限 b 为积分上限tol 为计算精度，缺省为0.001trace非0时，以动态点图的形式实现积分的整个过程。,-注意：调用quad函数时，先要建立一个描述被积函数f(x)的函数文件或语句函数。,例15.用自适应辛普森方法求积分,解：使用语
16、句函数（内联函数）,（高精度近似值为1.49364826）,clear;g=inline(exp(-x.2);,z=quad(g,-1,1),z=1.49364827,例16.用自适应辛普森方法求积分,（高精度近似值为1.49364826）,clear;g=inline(exp(-x.2);,z=quad(g,-1,1,1e-6),z=1.49364828,clear;g=inline(exp(-x.2);,z=quad(g,-1,1,1e-6,1),9-1.0000000000 5.43160000e-001 0.3198710950 11-1.0000000000 2.71580000e-
17、001 0.1290739989 13-0.7284200000 2.71580000e-001 0.1907965490 15-0.4568400000 9.13680000e-001 0.8538774475 17-0.4568400000 4.56840000e-001 0.4269534258 19-0.4568400000 2.28420000e-001 0.2024448313 21-0.2284200000 2.28420000e-001 0.2245087588 23 0.0000000000 4.56840000e-001 0.4269534258 25 0.00000000
18、00 2.28420000e-001 0.2245087588 27 0.2284200000 2.28420000e-001 0.2024448313 29 0.4568400000 5.43160000e-001 0.3198710950 31 0.4568400000 2.71580000e-001 0.1907965490 33 0.7284200000 2.71580000e-001 0.1290739989z=1.49364827606288,例17.用自适应辛普森方法求积分,解：建立函数文件,（高精度近似值为1.49364826）,function g=li12(x)g=exp(
19、-x.2);,z=quad(li12,-1,1),z=1.49364828,练习：用自适应辛普森方法法计算积分,6 二重积分,解：输入指令：,syms x y;int(int(x*y,y,2*x,x2+1),x,0,1),ans=1/12,syms x y;int(int(sin(pi*(x2+y2),y,-sqrt(1-x2),sqrt(1-x2),x,-1,1),解：积分区域可用不等式表成：,二重积分可化为二次积分,ans=sum(1/2*2(-3-2*_k1)*(-1)_k1*pi(2+2*_k1)/(1+_k1)*hypergeom(1/4,3/4,1/2,2+_k1,3/2+_k1,
20、-1/4*pi2)*2(1/2)/pi(1/2)*2(3/2+2*_k1)/(3/4+_k1)*gamma(5/2+2*_k1)/gamma(2+2*_k1)2+1/2*2(-2-2*_k1)*(-1)_k1*pi(2+2*_k1)*hypergeom(3/4,5/4,3/2,2+_k1,3/2+_k1,-1/4*pi2)*2(1/2)/pi(1/2)*2(2*_k1-1/2)/(1/4+_k1)/_k1*gamma(3/2+2*_k1)/gamma(2*_k1)/gamma(3+2*_k1),_k1=0.inf),syms r a;int(int(r*sin(pi*r2),r,0,1),a,
21、0,2*pi),输入命令：,结果中仍带有 int，表明 matlab 求不出这一积分，所以采用极坐标化为二次积分：,ans=2,7.数值计算矩形区域上二重积分,z=dblquad(f,a,b,c,d)求得二元函数f(x,y)的重积分，其中a,b为变量 x 的积分下、上限；c,d 为变量 y 的积分下、上限;z=tripquad(fun,a,b,c,d,e,f)求得三元函数fun(x,y,z)的重积分，其中a,b为变量 x 的积分下、上限；c,d 为变量 y 的积分下、上限;e,f 为变量 z 的积分下、上限,例20.计算重积分,解：,fun=inline(x.*exp(x.2+y.2),x,y
22、);dblquad(fun,0,2,-2,2),ans=,例21.计算重积分,解：,fun=inline(y.*sin(x)+z.*cos(x),x,y,z)triplequad(fun,0,pi,0,1,-1,1),ans=1.99999999436264,8 曲线积分,例22,解：,曲线的参数方程为,曲线积分可化为：,输入命令,ans=1/2,syms t;int(cos(t)*sin(t),0,pi/2),当曲线由参数方程给出,则曲线长为,做函数文件pmqxhc.m,function y=pmqxhc(x,y,t,a,b)y=int(sqrt(diff(x,t)2+diff(y,t)2),a,b);,例23,输入命令：,syms t a;x=a*(t-sin(t);y=a*(1-cos(t);s=pmqxhc(x,y,t,0,2*pi),pretty(simple(s)8 a,练习：,6.求区间等分数n=200,分别用梯形法和抛物线法编程，计算定积分,
