MATLAB第3章第3节.ppt
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 卫星轨道长度问题,椭圆的参数方程为 弧长的公式 椭圆长度,椭圆积分,无法用解析方法计算,讨论用数值方法来计算。,3.3.2 公式的导出,给定函数f(x)，关于积分，有Newton-Leibniz 公式 但是，在下列情况下，函数在离散点处给出；被积函数的原函数无法用初等函数表示；被积函数的原函数虽有初等函数表达式，但过于复杂；,就必须借助数值方法来求函数的积分,3.3.2 公式的导出,用数值方法近似地求一个函数 在区间(a,b)上的 定积分的基本思路，可以归结到定积分的定义,(3-4),3.3.2 公式的导出,取等距步长，当n充分大时，就是 的数值积 分，是第k小区间中x的取值，显然，取值不同，数值积分 的结果就不同。这种做法相当于用相对简单的阶梯函数(k=1,n)代替 作积分。实际上各种不同的数值积分方法就在于，研究用什么样的简单函数代替，使得既能保证结果有一定的精度，计算量又小。,3.3.2 公式的导出,牛顿-柯特斯(Newton-Cotes)公式 设 为给定的求积结点，将其作为插值结点，作 的拉格朗日插值多项式，然后，利用拉格朗日插值多项式替代 作积分，选取不同的多项式，就得到不同的求积公式。设，将 区间n等分，记，以这 个等距结点为求积结 点的插值型求积公式,通常被称为牛顿-柯特斯(Newton-Cotes)公式,3.3.2 公式的导出,常用的Newton-Cotes公式：n=0 分别用 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 辛普森公式的几何意义,3.3.2 公式的导出,复化的牛顿-柯特斯(Newton-Cotes)公式为构造高精度的数值积分公式，可以采用分段低次多项式替代整体高次多项式，这就导出了复化的牛顿-柯特斯公式，其基本思想是：先把积分区间分成一些长度较小的子区间，在每个子区间上使用低阶的牛顿-柯特斯公式，最常用的是下面的复化梯形公式和复化辛普森公式。设在(a,b)上，定积分表示曲线下的面积，我们先从图形上看看如何近似计算这块面积。,3.3.2 公式的导出,图3.7 定积分的复化矩形公式和复化梯形式,3.3.2 公式的导出,将(a,b)区间n等分，成为积分步长。记在每个小区间上用矩形面积近似 下面曲边梯形的面积，在整个区间(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段之和可得整个区间上的近似积分(3-8)3-8式称为复化辛普森求积公式（抛物线公式）。,3.3.3 误差估计和收敛性,误差估计 有了求积公式，如何度量它对原积分的近似程度呢？一种方式是考察 另 一种方式是用使得 的函数类的大小来度量。,人们称它为求积余项,引出了代数精确度的概念,3.3.3 误差估计和收敛性,（k阶代数精度）一个数值积分有k阶代数精度是指：如果当 是次数小于或等于k的多项式时，而对于 k+1次多项式，。即对任意次数不高于k次的 多项式，数值积分没有误差。,定义：,由拉格朗日插值多项式的性质可知，如上构造的求积公式的代数精确度至少是n.梯形公式：求积余项为(3-9)辛普森公式：求积余项为(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是由每个小区间中间点的函数值构成的向量如利用左矩形公式计算,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 应代入 此时步长相同且固定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.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.851506759375514e-004r2=6.851506759371073e-004r3=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+004l2=4.908996526863898e+004,
