《MATLAB第3章第3节.ppt》由会员分享，可在线阅读，更多相关《MATLAB第3章第3节.ppt（33页珍藏版）》请在三一文库上搜索。
1、3.3 数值积分,3.3.1 卫星轨道长度问题 3.3.2 公式的导出 3.3.3 误差估计和收敛性 3.3.4 用Matlab作数值积分 3.3.5 实验,3.3.1 卫星轨道长度问题,问题： 人造地球卫星轨道可视为平面上的椭圆。我国第一颗人造地球卫星近地点距地球表面439km ，远地点距地球表面2384km，地球半径为 6371km。 求: 该卫星的轨道长度。,3.3.1 卫星轨道长度问题,模型 a,b分别是长半轴和短半轴; 焦距为c， 地球半径为 r ， 近地点和远地点与地球表面的距离分别是 和 。,图3.4 卫星轨道的示意图,3.3.1 卫星轨道长度问题,椭圆的参数方程为 弧长的公式
2、椭圆长度,椭圆积分,无法用解析方法计算, 讨论用数值方法来计算。,3.3.2 公式的导出,给定函数f(x)，关于积分，有Newton-Leibniz 公式 但是，在下列情况下， 函数在离散点处给出； 被积函数的原函数无法用初等函数表示； 被积函数的原函数虽有初等函数表达式，但过于复杂；,就必须借助数值方法来求函数的积分,3.3.2 公式的导出,用数值方法近似地求一个函数 在区间(a,b)上的 定积分的基本思路，可以归结到定积分的定义,(3-4),3.3.2 公式的导出,取等距步长，当n充分大时， 就是 的数值积 分， 是第k小区间中x的取值，显然， 取值不同 ，数值积分 的结果就不同。 这种做
3、法相当于用相对简单的阶梯函数 (k=1,n) 代替 作积分。 实际上各种不同的数值积分方法就在于，研究用什么样的简单函数代替 ，使得既能保证结果有一定的精度，计算量又小。,3.3.2 公式的导出,牛顿-柯特斯(Newton-Cotes)公式 设 为给定的求积结点，将其作为插值结点，作 的拉格朗日插值多项式，然后 ，利用拉格朗日插值多项式替代 作积分，选取不同的多项式，就得到不同的求积公式。 设 ，将 区间n等分，记 ， 以这 个等距结点为求积结 点的插值型求积公式,通常被称为牛顿-柯特斯(Newton-Cotes)公式,3.3.2 公式的导出,常用的Newton-Cotes公式： n=0 分别
4、用 f (a)，f (b) 和 近似 f (x) 可得 几何解释：用以点 （左矩形公式）为顶点的矩形面积近似所要求的积分（矩形的面积）,右矩形公式,左矩形公式,中矩形公式,3.3.2 公式的导出,n=1 此时， ， 若用 f (a) 和 f (b) 的算术平均值近似 f ()，则可得 几何解释：用以点 为顶点的梯形面积近似所要求的积分（曲边梯形的面积）,梯形公式,图3.5 矩形公式、梯形公式的几何意义,3.3.2 公式的导出,3.3.2 公式的导出,n=2 此时， 求积公式为,通常称此公式为辛普森(Simpson)公式，也称 为抛物线公式。,3.3.2 公式的导出,图3.6 辛普森公式的几何意
5、义,3.3.2 公式的导出,复化的牛顿-柯特斯(Newton-Cotes)公式 为构造高精度的数值积分公式，可以采用分段低次多项式替代整体高次多项式，这就导出了复化的牛顿-柯特斯公式，其基本思想是：先把积分区间分成一些长度较小的子区间，在每个子区间上使用低阶的牛顿-柯特斯公式，最常用的是下面的复化梯形公式和复化辛普森公式。 设在 (a, b)上 ，定积分表示曲线下的面积，我们先从图形上看看如何近似计算这块面积。,3.3.2 公式的导出,图3.7 定积分的复化矩形公式和复化梯形式,3.3.2 公式的导出,将(a ,b)区间n等分， 成为积分步长。 记 在每个小区间上 用矩形面积近似 下面曲边梯形
6、的面积， 在整个区间(a ,b)内构成台阶形。,3.3.2 公式的导出,容易看出，两个台阶形面积分别为 在图3.7中，两个台阶形分别小于和大于所求面积故 3-5、 3-6就是计算定积分的复化矩形公式。,(3-5),(3-6),或将两者平均，则每个小区间上的小矩形 变为 小梯形，整个区间上的结果为,视 为结点，3-7式相当于用分段线性插值函数作为 的近似，称为复化梯形求积公式。,(3-7),3.3.2 公式的导出,3.3.2 公式的导出,用分段二次插值函数代替，记n=2m,k=0,1,m-1 在第k段的两个小区间上，用三个结点 作二次插值函数 ，然后积分，求m段之和可得 整个区间上的近似积分 (
7、3-8) 3-8式称为复化辛普森求积公式（抛物线公式）。,3.3.3 误差估计和收敛性,误差估计 有了求积公式，如何度量它对原积分的近似 程度呢？ 一种方式是考察 另 一种方式是用使得 的函数类的 大小来度量。,人们称它为求积余项,引出了代数精确度的概念,3.3.3 误差估计和收敛性,（k阶代数精度） 一个数值积分有k阶代数精度是指：如果当 是次数小于或等于k的多项式时， ，而对于 k+1次多项式， 。即对任意次数不高于k次的 多项式 ，数值积分没有误差。,定义：,由拉格朗日插值多项式的性质可知，如上构 造的求积公式的代数精确度至少是n. 梯形公式：求积余项为 (3-9) 辛普森公式：求积余项
8、为 (3-10),1阶代数精度,3阶代数精度,3.3.3 误差估计和收敛性,收敛性 若对 的某个数值积分 有 （非零常数），则称 是 p阶收敛的。 按此定义可以判断梯形公式是2阶收敛的，类似地，辛普森公式是4阶收敛的。,定义：,3.3.3 误差估计和收敛性,（p阶收敛的）,3.3.4 用Matlab作数值积分,对于向量x，cumsum(x)返回一结果向量，此向量x的第n个元素为向量的前个元素之和，如 x=1 1 1 1; I=cumsum(x) I = 1 2 3 4,矩形求积指令,cumsum(x),x是由每个小区间左端点的函数值构成的向量 x是由每个小区间右端点的函数值构的向量 x是由每个
9、小区间中间点的函数值构成的向量 如利用左矩形公式计算,cumsum(x)*h,其中h为子区间步长,左矩形公式,右矩形公式,中矩形公式,3.3.4 用Matlab作数值积分,计算矩形积分公式,3.3.4 用Matlab作数值积分, x=linspace(0,pi,50); h=1:49； t(h)=x(h); y=sin(t); T=cumsum(y)*pi/(49); I=T(49) I = 1.9993,（指令x=linspace(a,b,n): 在a,b区间中的n个等分点（包括端点）构成的向量）,梯形求积指令trapz,用梯形方法计算Y的积分近似值。对于向量Y，Y 应代入 此时步长相同且固
10、定h=1 。 使用梯形法计算Y对X的积分。其中输入量 ：,；,Z=trapz(Y),Z=trapz(X,Y),3.3.4 用Matlab作数值积分,3.3.4 用Matlab作数值积分,如 x=linspace(0,pi,50); format long y=sin(x);T=trapz(x,y) T = 1.99931484932406,梯形求积指令trapz 计算,辛普森公式求积指令quad,用simpson公式求函数fun(x)在a,b的积分近似值,自动选取步长,相对误差为,输出积分值。 同上,但相对误差为tol。,quad(fun,a,b),quad(fun,a,b,tol),3.3.
11、4 用Matlab作数值积分,几种求积指令的比较,计算积分 （很显然这个积分的值是2） 编写M文件如下： x=linspace(0,pi,50); format long y=sin(x); T1=cumsum(y(1:49)*pi/(49); %左矩形法 T1=T1(49) T2=trapz(x,y); %梯形法,3.3.4 用Matlab作数值积分,3.3.4 用Matlab作数值积分,T3=quad(sin,0,pi); %辛普森法 % r1 r2 r3是几种方法的绝对误差 r1=abs(T1-2); r2=abs(T2-2); r3=abs(T3-2);,计算结果： r1 = 6.85
12、1506759375514e-004 r2= 6.851506759371073e-004 r3= 3.601569042999131e-009,3.3.5 实验,利用梯形公式和辛普森公式计算3.3.1中 卫星轨道长度问题，编写M文件如下： function y=x5(t) a=8755;b=6810; y=sqrt(a2*sin(t).2+b2*cos(t).2); 在命令窗口键入以下程序： t=0:pi/10:pi/2; y1=x5(t); l1=4*trapz(t,y1);,定义一个函数,调用这个函数,3.3.5 实验,l2=4*quad(x5,0,pi/2,1e-6); format long,用长格式输出结果： l1=4.908996526785276e+004 l2= 4.908996526863898e+004,