5.4节_数字积分和微分方程数值解

5.4节数值积分和微分方程数值解一．数值定积分求面积•【例5-4-1】用数值积分法求由，y=0,x=0与x=10围成的图形面积，并讨论步长和积分方法对精度的影响。•解:◆原理用矩形法和梯形法分别求数值积分并作比较，步长的变化用循环语句实现。MATLAB中的定积分有专门的函数QUAD,QUADL等实现。为了弄清原理，我们先用直接编程的方法来计算，然后再介绍定积分函数及其调用方法。设x向量的长度取为n，即将积分区间分为n-1段，各段长度为。算出各点的，则矩形法数值积分公式为：2115yx(1,2,,1)ixin(1,2,,1)iyin11niiisyx矩形和梯形定积分公式•梯形法的公式为：•比较两个公式，它们之间的差别只是。•在MATLAB中，把向量中各元素叠加的命令是sum。把向量中各元素按梯形法叠加的命令是trapz。梯形法的几何意义是把被积分的函数的各计算点以直线相联，形成许多窄长梯形条，然后叠加，我们把两种算法都编入同一个程序进行比较。11111222nniiniiiiiyyyyqxyx10.5nyy求面积的数值积分程序exn541fordx=[2,1,0.5,0.1]%设不同步长x=0:.1:10;y=-x.*x+115;%取较密的函数样本plot(x,y),holdon%画出被积曲线并保持x1=0:dx:10;y1=-x1.*x1+115;%求取样点上的y1%用矩形（欧拉）法求积分，注意末尾去掉一个点n=length(x1);s=sum(y1(1:n-1))*dx;q=trapz(y1)*dx;%用梯形法求积分stairs(x1,y1),%画出欧拉法的积分区域plot(x1,y1)%画出梯形法的积分区域[dx,s,q],pause(1),holdoff，end程序exn541运行结果程序运行的结果如下：步长dx矩形法解s梯形法解q29108101865815.5841.25816.25.1821.65816.65•用解析法求出的精确解为2450/3=816.6666...。•dx=2时矩形法和梯形法的积分面积见图5-4-1.。在曲线的切线斜率为负的情况下，矩形法的积分结果一定偏大，梯形法是由各采样点联线包围的面积，在曲线曲率为负（上凸）时，其积分结果一定偏小，因此精确解在这两者之间。由这结果也能看出，在步长相同时，梯形法的精度比矩形法高。矩形法数字积分的演示程序rsums•MATLAB中有一个矩形法数字积分的演示程序rsums，可以作一个对比。键入rsums('115-x.^2',0,10)•就得到右图。图中表示了被积函数的曲线和被步长分割的小区间，并按各区间中点的函数值构成了各个窄矩形面积。用鼠标拖动图下方的滑尺可以改变步长的值，图的上方显示的是这些矩形面积叠加的结果。0246810020406080100115-x.^2:816.99218816MATLAB内的数值定积分函数•在实际工作中，用MATLAB中的定积分求面积的函数quad和quadl可以得到比自编程序更高的精度，因为quad函数用的是辛普生法，即把被积函数用二次曲线逼近的算法，而quadl函数采用了更高阶的逼近方法。它们的调用格式如下：Q=QUADL(FUN,A,B,TOL)其中，FUN是表示被积函数的字符串，A是积分下限，B是积分上限。TOL是规定计算的容差，其默认值为1e-6•例如，键入S=quad('-x.*x+115',0,10)得到S=8.166666666666666e+002二．求两条曲线所围图形的面积【例5-4-2】。设计算区间[0,4]上两曲线所围面积。•解：◆原理：先画出图形，•dx=input('dx=');x=0:dx:4;•f=exp(-(x-2).^2.*cos(pi*x));•g=4*cos(x-2);•plot(x,f,x,g,':r')2(2)cos(),()4cos(2)xxfxegxx得到右图。从图上看到，其中既有f(x)g(x)的区域，也有f(x)g(x)的区域，求两条曲线所围图形的面积(1)若要求两曲线所围总面积（不管正负），则可加一条语句s=trapz(abs(f-g))*dx,在dx=0.001时，得到s=6.47743996919702若要求两曲线所围的f(x)g(x)的正面积，则需要一定的技巧.◆方法一。先求出交点x1，再规定积分上下限。x1=fzero('exp(-(x-2).^2.*cos(pi*x))-4*cos(x-2)',1)%把积分限设定为0~x1，求出积分结果再乘以2：x=0:dx:x1;f=exp(-(x-2).^2.*cos(pi*x));g=4*cos(x-2);s1=2*trapz(abs(f-g))*dx在设定dx=0.001时，得到s1=2.30330486000857求两条曲线所围图形的面积(2)方法二。调用MATLAB中求面积函数quad。这里的关键是建立一个函数文件，把e1=f(x)-g(x)0的部分取出来。•利用逻辑算式(e10)，它在e10处取值为1，在e10处则为零。让逻辑函数(e10)与e1作元素群乘法，正的e1将全部保留，而负的e1就全部为零。因此编出子程序exn542f.m如下：•functione=exn542f(x)•e1=exp(-(x-2).^2.*cos(pi*x))-4*cos(x-2);•e=(e1=0).*e1;•将它存入工作目录下。于是求此积分的主程序语句为：•s2=quad('exn542f',0,4)•得到的结果为：•s2=2.30330244952618三．求曲线长度【例5-4-3】设曲线方程及定义域为：用计算机做如下工作：(a)按给定区间画出曲线，再按n=2,4,8份分割并画出割线。(b)求这些线段长度之和，作为弧长的近似值。(c)用积分来估算弧长，并与用割线计算的结果比较。解：◆原理：先按分区间算割线长度的方法编程，然后令分段数不断增加求得其精密的结果，最后可以与解析结果进行比较。因此编程应该具有普遍性，能由用户设定段数，并在任何分段数下算出结果。2()1,11fxxx求曲线长度的程序exn543n=input('分段数目n='),%输入分段数目x=linspace(-1,1,n+1);%设定x向量y=sqrt(1-x.^2);%求y向量plot(x,y),holdon%绘图并保持Dx=diff(x);%求各段割线的x方向长度%x向量长度为n+1，Dx是相邻x元素的差，其元素数为nDy=diff(y);%Dy是相邻两个y元素的差Ln=sqrt(Dx.^2+Dy.^2);%求各割线长度L=sum(Ln)%求n段割线的总长度程序exn543的运行结果•程序运行后得到图5-32，在不同的n下，其数值结果为：•n=2,L=2.82842712474619•n=4,L=3.03527618041008•n=8,L=3.10449606777484•n=1000L=3.14156635621648•我们已经可以大致猜测出它将趋向于π，精确的极限值可用下列符号数学语句导出。•x=symsx,y=sqrt(1-x^2),L=int(sqrt(1+diff(y)^2),-1,1)•这个程序其实有相当的通用性，不同的被积函数，只要改变其中的一条函数赋值语句，并相应地改变自变量的赋值范围就行了。四．求旋转体体积【例5-4-4】求曲线与x轴所围成的图形分别绕x轴和y轴旋转所成的旋转体的体积。解：原理：由于旋转对称性，在圆周方向的计算只要乘以圆周长度，不需要积分运算。因此旋转体的体积计算实际上就退化单变量求积分。程序如下：%先画出平面图形dx=input('dx=');x=0:dx:pi;g=x.*sin(x).^2;plot(x,g)求筒形旋转体体积(a)。绕y轴体积用薄圆柱筒形体作为微分体积单元，其半径为x，厚度为dx，高度为g(x)，其立体图见图5-34左，此筒形单元的截面积为g*dx，薄环的微体积为：dv=2*pi*x*dx*g，旋转体的体积为微分体积单元沿x方向的和，键入：v=trapz(2*pi*x.*g*dx)得：v=27.53489480036561求盘形旋转体体积(b)。绕x轴体积它绕x轴旋转形成一个薄圆盘，其厚度为dx，而半径为g(x)。所以此薄盘体的微体积为：dv1=pi*g.^2.*dx，旋转体的体积为微分体积单元沿x方向的和：v1=trapz(pi*g.^2*dx)得：v1=9.86294784774497用符号数学工具箱的程序•精确的理论结果可用符号数学工具箱函数求得如下：•symsx,g=x*sin(x)^2;•v1t=int(pi*g^2,0,pi),double(v1t)•v1t=1/8*pi^4-15/64*pi^2•ans=9.86294784774499•大多数的定积分并不会有理论的解析结果，所以这样的验证一般是不必要的。五．多重积分【例5-4-5】计算二重积分积分区域Ω为由x=1,y=x及y=0所围成的闭合区域.22xydxdy•解:◆原理先画出积分区域，在任意x处取出沿y向的一个单元条，其宽度为dx，而高度为y=x，所以y是一个数组。其上的被积函数f也是一个数组，沿y向的积分可用trapz(f)完成，得到s1(k)，它是随x而变的。用for循环求出所有的s1(k)。再沿x方向用trapz函数积分。MATLAB的数组运算可以代替一个for循环，所以二重积分只用了一组for语句。二重积分的MATLAB程序exn545clear,formatcompactfill([0,1,1,0],[0,0,1,0],'y'),hold%画出积分区域fill([0.55,0.6,0.6,0.55,0.55],[0,0,0.6,0.55,0],'r')%画出单元条dx=input(‘步长dx=‘);dy=dx;x=0:dx:1;lx=length(x);fork=1:lxx1=(k-1)*dx;y1=0:dy:x1;f=x1.^2+y1.^2;s1(k)=trapz(f)*dy;ends=trapz(s1)*dx用MATLAB函数求二重积分(1)运行的数值结果在步长dx=0.01时为:s=0.3334另一种方法是利用MATLAB中现成的二重积分函数dblquad，其调用格式为：Q=DBLQUAD(FUN,XMIN,XMAX,YMIN,YMAX,TOL)•其中FUN是x,y的函数，接下来的四个变元是四个积分限，其中前两个对应于x，后两个对应于y，TOL为允许误差（默认值为1.e-16）。这四个积分限只许用常数代入，可见dblquad函数只能用于积分区域为矩形的情况。•解决的方法之一是仍用矩形区积分，但把不属于积分区域内的函数置成零，其方法与上题有些类似。用MATLAB函数求二重积分(2)•在图示的积分区域中，对角线左上方的白色区域满足y-x0，逻辑式(y-x0)在此区域均等于零,而在灰色区域内为一。将它与被积函数作元素群相乘，就构成了一个新的被积函数，它与原来函数的差别是把灰色积分区外的函数置为零。这样就可以按矩形区域调用dblquad函数了。键入：Q=dblquad('(x.^2+y.^2).*(y-x0)',0,1,0,1)•得到Q=0.33332245532028三重积分的计算【例5-4-6】计算三重积分积分区域Ω为由x=1,y=x,z=xy及三个坐标面所围成的区域.解:◆方法先画出积分区域图5-36,这个区域在xy平面上的投影与图5-35相仿，只是增加了z方向的高度，从而构成了一个三维的实体。23xyzdxdydz先画出积分区域。这个区域在xy平面上的投影上例相仿，只是增加了z方向的高度，从而构成了一个三维的实体。程序exn546a用来画这个立体空间。•x=1,y=x都是沿z向的柱面，本题还