1、数学实验3,插值与数值积分(2),本次课的主要内容,数值积分的四种方法1、矩形公式sum 或 cumsum2、梯形公式trapz3、辛普森公式quad4、高斯公式编程2：gaussinteg Gauss-Lobatto方法 quadl,编程1：simp,3.3 数值积分,数值积分:用数值的方法近似地求一个定积分,(教材P56,参考书P222),数值分析中主要介绍三种等距节点的求积公式(牛顿-科茨公式) :,1、矩形公式（k=0）,h,h=(b-a)/n,fk=f(xk),2、梯形公式（k=1）,3、辛甫森公式（k=2）,已知n+1对节点数据(xi , yi)(i=0,1,n)，求积分.,tra
2、pz(y) 按梯形公式计算定积分（单位步长）。 trapz(x,y) x , y同长度，输出 y 对 x 的按梯形公式计算的积分 （变步长）。,quad(fun,a,b) 用辛甫森（2阶）公式计算函数fun在区间 a, b的积分，自动选择步长。quad(fun,a,b,tol) 与上同，但指定了相对误差 tol。,1 MATLAB作数值积分,一 用于数值积分的几种命令：,sum(x) 输入数组x，输出为x的和，用于矩形公式求积分。 通常还要乘以等分小区间的长度(b-a)/ncumsum(x) 返回与x一样长的向量。此向量的第n个元素为x的 前n个元素之和.用于矩形公式求积分，同样要乘 (b-a
3、)/n（参考书P223）,quadl(fun,a,b,tol) 用自适应Gauss-Lobatto公式计算，精度更高。,function s=simp(x,y)if mod(length(x),2)=0 error(数据点必须为奇数)endn=length(x);m=(n-1)/2;h=(x(n)-x(1)/2/m;s1=0; s2=0;for i=1:m s1=s1+y(2*i);endfor j=1:m-1 s2=s2+y(2*j+1);ends=(y(1)+y(n)+4*s1+2*s2)*h/3,根据辛甫森求积公式,当被积函数不是解析表示时,比如离散数据表表示的函数通常就用这个函数按辛甫
4、森公式计算积分。,将a,b区间等分成2m个小区间：,每2个区间构成一个组,即：,每一组有3个节点，如第一组：,由3个节点构造一个拉格朗日2次多项式：,返回,令f(x)= xk ，用(11)式计算,我们不妨只考虑,二 高斯（Gauss)求积公式,各种近似求积公式都可以表示为,(11),若对于k=0,1,.,m 都,有In = I ，而当k =m+1时，In I ，则称In 的代数精度为m 。梯形公式代数精度为 1，辛甫森公式的代数精度为 3。,下面介绍的是取消对区间等分的限制，n 给定后同时确定节点xi和系数Ai，使代数精度尽可能高的所谓高斯公式。,我们先考虑节点数为2 而使用(11)计算的积分
5、近似值有代数精度为3。,而构造代数精度为3的形如,G2=A1 f(x1)+ A2 f(x2),(12),的求积公式。,成立，依次将f(x)= 1, x, x 2, x 3代入，即可得到确定A1,A2 ,x1 ,x2的方程组。,如何选择节点xi 和系数Ai ，使(11)计算的精度更高？,对于 f (x)=1, x, x 2, x 3，应该有,将 f (x) =1, x , x 2 , x 3 依次代入，得：,提高精度可以通过增加节点数n ，(11)的代数精度可达到2n-1,但增加了解高维线性方程组的难度，实用价值不大；另一方法是将区间分小，在小区间上用G2。,n=2的高斯公式为：,将区间(a,
6、b)作m等分，记h=(b-a)/m，xk=a+kh, k=0,1,.,m, 作变换,其中,常用的高斯公式就是：,下面我们来编写M文件，应用高斯公式计算定积分。,将x xk-1, xk化为t-1, 1,function s=gaussinteg(f,a,b,m)s=0;h=(b-a)/m;k=1:m+1;x(k)=a+(k-1)*h;j=1:m;z1(j)=(x(j)+x(j+1)/2+h/2/sqrt(3);z2(j)=(x(j)+x(j+1)/2-h/2/sqrt(3);for j=1:m s=s+feval(f,z1(j)+feval(f,z2(j); ends=s*h/2;,functi
7、on y=jifenfun(x)y=exp(x).*sin(x);return,定义高斯积分函数M文件：,s=gaussinteg(jifenfun, . 0,2*pi,1000)s = -267.2458,在命令窗口运行,定义任何一个被积函数,根据,Gauss-lobatto是改进的高斯积分方法,采取自适应求积方法,三 用MATLAB实现定积分计算：, 矩形公式与梯形公式,z1 = 0.9602z2 = 1.0388z11 = 0.9602z12 = 1.0388z3 = 0.9995,h=pi/40;x=0:h:pi/2;y=sin(x);z1=sum(y(1:20)*h,z2=sum(y
8、(2:21)*h,z=cumsum(y)*h, z11=z(20)*h,z12=(z(21)-z(1)*h,z3=trapz(y)*h,z=0 0.0062 0.38020.0184 0.44380.0368 0.51070.0611 0.58070.0911 0.65330.1268 0.72800.1678 0.80430.2140 0.88190.2650 0.96020.3205 1.0388, 用辛甫森公式求定积分,z4=quad(sin,0,pi/2)z5=quadl(sin,0,pi/2)z4 = 1.0000z5 = 1.0000, 高斯积分公式求定积分,s=gaussinte
9、g(sin, 0, pi/2,1000)s = 1.0000,前面已经保存了文件gaussinteg.m在搜索路径中，即可直接调用。这里被积函数为内部函数，无需另外定义。,一 求卫星轨道长度,我国第一颗人造地球卫星轨道为平面椭圆曲线，近地点距地球表面439km，远地点距地面2384km，地球半径为6371km，求卫星轨道的长度。,椭圆的参数方程为：,椭圆长度由以下椭圆积分表示：,a=7782.5 kmb=7721.5 km,2 数值积分应用问题举例,我们用梯形公式和辛甫生公式来求卫星轨道的长度。,t=0:pi/10:pi/2;y1=track1(t);format long el1=4*tra
10、pz(t,y1)l2=4*quad(track1,0,pi/2,1e-6),function y=track1(t)a=7782.5;b=7721.5;y=sqrt(a2*sin(t).2+. b2*cos(t).2);,l1 = 4.908996526785276e+004l2 = 4.908996531830460e+004,先创建一个函数M文件：track1.m,再建立一个脚本文件,tracklong.m,运行结果：,广义积分、二维数值积分,1、无穷区间截断误差：blquadtriplequad,p65 3 , 4 2) ,11,3、要求用一个脚本文件把三种公式计算的结果和计算出的精确结
11、果用矩阵表示，在输出中能直观比较。4、用梯形公式计算要求用不同的步长计算作比较，用辛甫森公式计算要求用不同的精度计算 作比较。设计好M文件保存输出结果便于比较。,请将simp.m, gaussinteg.m, 编辑好保存在自己的文件夹中，以便随时调用。,上机作业,四 蒙特卡罗方法(用随机模拟方法计算数值积分),例、向图中边长为 1 的正方形 里随机投 n 块石头，落在四分之一圆内的石头为 k 块，根据几何概率，则四分之一圆的面积的近似值就是k / n，即,应用：先由计算机随机产生n个点的坐标（xi , yi )(i=1,2,n),其中xi ,yi 分别为a，b和c，d区间上的均匀分布随机数，然
12、后记录n个点中落在区域S(满足 )中的数目k.,x=a+rand(1,n)*(b-a);y=c+rand(1,n)*(d-c);k=0;for i=1:n if y(i)=f(x(i) k=k+1; endends=k/n*(b-a)*(d-c),1、随机投点法：,若随机变量X 的概率分布密度是 p(x), a x b，则,特别当X为a, b区间均匀分布的随机变量时，p(x)=1/(b-a), (axb),这两种用随机模拟的方式求积分近似值的方法,注2、无须作随机数 yi f (xi ) 比较。,注1、xk(k=1,2,n)为a,b 的随机数,2 均值估计法：,蒙特卡罗方法,通常还要先设法使被
13、积函数在0, 1范围取值再应用随机投点法。,function z=motc2(f,a,b,n)x=rand(1,n);y=rand(1,n); k=0; for i=1:n if y(i)z=motc2(sin,0,pi/2,10000),function z=motc1(f,a,b,n)l=b-a; z=0;for k=1:n t=rand(1); z=z+feval(f,a+l*t);endz=z*(b-a)/n;,z=motc1(sin,0,pi/2,1e4),z = 1.00130915270174,z = 0.98897336735007,3、蒙特卡罗方法的通用函数与调用格式,均值估
14、计法,随机投点法,(设0 f(x) 1),4、运用蒙特卡罗方法计算重积分,a、随机投点法,function s=motc3(f,x,y,z,n)If nargin5 n=10000 ; end;k=0;xh=x(2)-x(1);yh=y(2)-y(1);zh=z(2)-z(1);for i=1:n a=x(1)+xh*rand; b=y(1)+yh*rand; c=z(1)+zh*rand; u=feval(f,a,b,c); if u=0 k=k+1;end;ends=k/n*xh*yh*zh;,例、求半径为 的半球的体积。,计算二重积分：,设xi, yi(i =1,n)是相互独立，分别为a
15、, b和c,d的随机数，判断每个点(xi, yi) 是否落在W域内，将落在 W域内的m个点记作(xk, yk) ,k=1,m，则,蒙特卡罗方法:优点:计算简单,可以计算重积分;缺点:计算量大，精度较差，结果随机,b、均值估计法,3 数值微分,一 数值微分公式 当函数以离散数值形式给出时，用离散方法近似地计算函数y= f (x)在某点的导数值，这就是数值微分。,导数有以下三个近似公式：,它们依次称为前差公式、后差公式和中点公式。这些公式都明确的几何意义：,y=f(x),B,A,C,T,a-h a a+h,x,y,o,以上三个公式分别表示了线段AB、CA和CB的斜率。前、后差公式误差为O(h)中点
16、公式误差O(h2),后两个公式是由二次插值得到，目的是保持在端点处的精度O(h2),高阶导数的近似公式一般要用插值多项式得到，下面是二阶导数的近似公式：,diff(x): x2-x1 , x3-x2, xn-xn-1,二 应用实例人口增长率,已知20世纪美国人口统计数据如下，试计算表2 中这些年份的人口增长率。,又已知某地区20世纪70年代的人口增长率如表3，且1970年人口为210（百万），试估计1980年的人口。,1 记时刻t 的人口为x(t)，人口相对增长率为,，记,19001990年的人口依次为 xk ，(k=0,1,.,n)，年增长率为rk 。,由三点公式可以得到：,x=76.0 9
17、2.0 106.5 123.2 131.7. 150.7 179.3 204.0 226.5 251.4;r(1)=(-3*x(1)+4*x(2)-x(3)/20/x(1);for k=2:9 r(k)=(x(k+1)-x(k-1)/20/x(k);endr(10)=(x(8)-4*x(9)+3*x(10)/20/x(10);r,r =0.0220 0.0166 0.0146 0.0102 0.0104 0.0158 0.01490.0116 0.01050.0104,2 某地区20世纪70年代人口增长率如下表：,x0=210;r=0.87 0.85 0.89 0.91 0.95 1.10/100;y=trapz(r)*2; %步长为2x=x0*exp(y),x = 230.1676,p93 5.2 ) b，d 5) 10),10) 概率密度函数是,其中x, y分别为x方向和y方向的均方差,为偏差在x和y之间的相关系数以100m为一个单位，目标区x2+y2=1采用蒙特卡罗方法之均值估计法，参考P91的4.2射击命中概率,11)用蒙特卡罗法计算重积分。(思考),实验内容与要求,
