实验6 数值积分.doc

实验六 数值积分(Quadrature)6.1 实验目的掌握数值积分的代数精度与收敛阶的含义，会利用matlab求解符号积分和数值积分，并会借助数学软件Matlab求解一些简单的实际问题。6.2 实验内容1、Matlab中求解符号积分和数值积分的方法； 2、建立飞船的轨道周长和手的面积等实际问题的数学模型，并借助数学软件Matlab求解6.3 实验步骤6.3.1 函数表达式已知时的积分Matlab中求积分的方法包括符号积分（精确计算）和数值积分两种。前者只有一个命令int（注：求导数的命令为diff）,后者命令较多，详见表1。表1 Matlab中数值积分函数一览表函数名调用格式方法说明适应范围（详见Help）quadq=quad(fun,a,b,tol)自适应Simpson 法被积函数不光滑，低精度要求quadlq=quadl(fun,a,b,tol)自适应 Lobatto 法被积函数光滑，高精度要求 quadgkq,errbnd=quadgk(fun,a,b,tol)自适应Gauss- Kronrod 法震荡型被积函数，高精度要求，积分限可为无穷trapzq=trapz(X,Y)梯形法离散数据点积分dblquadq = dblquad(fun,xmin,xmax,ymin,ymax)调用一维数值积分方法矩形区域上的二重积分quad2dq = quad2d(fun,a,b,c,d)平面区域上的二重积分triplequadtriplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax)调用一维数值积分方法三维长方形区域上的三重积分教师示范1：符号积分int是matlab中符号工具箱的一个命令，它能精确计算不定积分和定积分，微积分教材中的积分一般均可用该种方法积出。例1 求不定积分。Matlab命令序列为：syms xs=int(-2*x/(1+x2)2)运行结果为：1/(x2 + 1)，即。验算所得结果，即验证，Matlab命令为：diff(s)，运行结果为：-(2*x)/(x2 + 1)2.其余见文件：integrate_symbolic.m。教师示范2：数值积分例2 求不定积分Matlab命令为：format longq1 = quad(x) exp(-x.2/2),0,2)q2 = quadl(x) exp(-x.2/2),0,2)q3 = quadgk(x) exp(-x.2/2),0,2)结果分别为1.196288044999181，1.196288103140270和1.196288013322608。实际上，由此，函数quadgk比quad、quadl的计算精度要高。6.3.2 离散数据的数值积分当函数的表达式未知，但已知函数在一系列节点处的函数值时，的计算有三种方法。方法1，利用表1中的梯形法；方法2，先用样条插值函数逼近未知函数，然后利用表1中的命令求该样条函数的积分；方法3，如果所求的图形可视为以点为顶点的多边形，则可利用多边形的面积公式（1） （1）直接求解。具体见例3。例3 用不同方法计算不定积分，比较不同方法的精度。程序如下。% 比较不同方法求定积分的精度a = 0; b = 10;syms x;exact = vpa( int(cos(x),a,b), 20); % 符号定积分，小数点20位的精确值 q = zeros(7,1);q(1) = quad(cos,a,b);q(2) = quadl(cos,a,b);q(3) = quadgk(cos,a,b);% 离散数据后，用梯形法求积分t = a:0.1:b; y = cos(t); q(4) = trapz(t, y); % 用三次样条近似函数后，求积分pp = spline( t, y );q(5) = quad(t)ppval(pp,t),a,b);q(6) = quadl(t)ppval(pp,t),a,b);q(7) = quadgk(t)ppval(pp,t),a,b); %disp(精确值)disp(exact)format long;disp(近似值)disp(q);disp(各种方法的误差)for i = 1:7 fprintf(1,%12.10en,double( q(i)-exact ) );end以上程序运行结果见表2。表2 Matlab中不同数值积分方法精度比较函数结果方法说明误差int-0.5440211108893698134符号积分，精确解quad-0.544021121312844自适应Simpson 法-1.04e-008quadl-0.544021110988547自适应 Lobatto 法-9.92e-011quadgk-0.544021110889370自适应Gauss- Kronrod 法-1.83e-016trapz-0.543567684387146梯形法4.53e-004splinequad-0.544021095007636三次样条近似1.59e-008quadl-0.5440209206203261.91e-007quadgk-0.5440209978193161.13e-007详见文件：integrate_numeric.m。6.3.2 学生练习学生练习1：网球拍的拍面外形可近似看成椭圆，如图1所示。试计算该网球拍拍面外形的周长。首先给出求解椭圆周长的积分公式，然后数值求解。图1. 网球拍的拍面据中国载人航天工程网报道1，中国第一艘载人航天飞船神舟五号的初始轨道为近地点200公里、远地点350公里的椭圆轨道，试计算该椭圆轨道的长度。设地球为半径6371公里的球形。（选做）你可以通过查阅权威网站(如中国航天网)提供的数据，计算嫦娥二号绕月飞行的椭圆轨道的长度。学生练习2：实验5插值法的学生实验练习2中已绘制出你的手，试计算你的手的面积。图2. 一只手根据文献2，“用病人的手估计估计烧伤面积已成为诊断烧伤面积的公认方法之一”，并指出“男性手掌的面积大致占人体表面积的0.81%，女性则大致为0.67%”。试由此估计你的体表面积。提示：手的面积的求法有多种。一种最简单的方法是利用多边形的面积。可选的方法是定积分或二重积分，请同学们认真思考。给出不少于两种的计算方法。6.4 思考题（选做）1用数值积分方法计算含参变量的积分，并将你的计算结果与Matlab中的内置函数beta(z,w)和gamma(x)的结果进行比较。前者可参见integrate_numeric.m。2查阅资料学习如何用Monte-Carlo算法计算二重积分。然后，分别用不同的，计算积分，其中，对于每