第十讲 数值积分和微分方程数值解.ppt

第十讲 数值积分和微分方程数值解 一．数值定积分求面积 【例5-4-1】 用数值积分法求由 ，y=0, x=0与x=10围成的图形面积，并讨论步长和积分方法对精度的影响。 解: ◆原理 用矩形法和梯形法分别求数值积分并作比较，步长的变化用循环语句实现。MATLAB中的定积分有专门的函数QUAD,QUADL等实现。为了弄清原理，我们先用直接编程的方法来计算，然后再介绍定积分函数及其调用方法。设x向量的长度取为n，即将积分区间分为n-1段， 各段长度为 。算出各点的 ，则矩形法数值积分公式为： 矩形和梯形定积分公式 梯形法的公式为： 比较两个公式，它们之间的差别只是 。 在MATLAB中，把向量中各元素叠加的命令是sum。把向量中各元素按梯形法叠加的命令是trapz。梯形法的几何意义是把被积分的函数的各计算点以直线相联，形成许多窄长梯形条，然后叠加，我们把两种算法都编入同一个程序进行比较。 求面积的数值积分程序exn541 for dx=[2,1,0.5,0.1] % 设不同步长 x=0:.1:10;y=-x.*x+115; % 取较密的函数样本 plot(x,y),hold on % 画出被积曲线并保持 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); % 画出欧拉法的积分区域 hold on; plot(x1,y1) ; % 画出梯形法的积分区域 [dx,s,q],pause(1), hold off; end 程序exn541运行结果 程序运行的结果如下： 步长dx 矩形法解s 梯形法解q 2 910 810 1 865 815 .5 841.25 816.25 .1 821.65 816.65 用解析法求出的精确解为2450/3=816.6666...。 dx=2时矩形法和梯形法的积分面积见图5-4-1.。在曲线的切线斜率为负的情况下，矩形法的积分结果一定偏大，梯形法是由各采样点联线包围的面积，在曲线曲率为负（上凸）时，其积分结果一定偏小，因此精确解在这两者之间。由这结果也能看出，在步长相同时，梯形法的精度比矩形法高。 矩形法数字积分的演示程序rsums MATLAB中有一个矩形法数字积分的演示程序rsums，可以作一个对比。键入 rsums(115-x.^2,0,10) 就得到右图。图中表示了被积函数的曲线和被步长分割的小区间，并按各区间中点的函数值构成了各个窄矩形面积。用鼠标拖动图下方的滑尺可以改变步长的值，图的上方显示的是这些矩形面积叠加的结果。 MATLAB内的数值定积分函数 在实际工作中，用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) 求两条曲线所围图形的面积(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))