数学实验3 插值与数值积分(2)

本次课的主要内容 编程1：simp 数值积分的四种方法 1、矩形公式 sum 或 cumsum 2、梯形公式 trapz 3、辛普森公式 quad 4、高斯公式 编程2：gaussinteg Gauss-Lobatto方法 quadl 编程1：simp

3.3 数值积分 (教材P54,参考书P222) 数值积分:用数值的方法近似地求一个定积分 已知n+1对节点数据(xi , yi)(i=0,1,…n)，求积分. 数值分析中主要介绍三种等距节点的求积公式(牛顿-科茨公式) : x1 x2 x3 xi xi+1 xn-1 f1 f2 f3 fi fi+1 f0 fn fn-1 x y o 1、矩形公式（k=0） a b h=(b-a)/n fk=f(xk) h 2、梯形公式（k=1） 3、辛甫森公式（k=2）

§1 MATLAB作数值积分 一 用于数值积分的几种命令： sum(x) 输入数组x，输出为x的和，用于矩形公式求积分。 一 用于数值积分的几种命令： sum(x) 输入数组x，输出为x的和，用于矩形公式求积分。 通常还要乘以等分小区间的长度(b-a)/n cumsum(x) 返回与x一样长的向量。此向量的第n个元素为x的 前n个元素之和.用于矩形公式求积分，同样要乘 (b-a)/n（参考书P223） quad('fun',a,b) 用辛甫森（2阶）公式计算函数fun在区间 [ a, b]的积分，自动选择步长。 quad('fun',a,b,tol) 与上同，但指定了相对误差 tol。 trapz(y) 按梯形公式计算定积分（单位步长）。 trapz(x,y) x , y同长度，输出 y 对 x 的按梯形公式计算的积分 （变步长）。 quadl(‘fun’,a,b,tol) 用自适应Gauss-Lobatto公式计算，精度更高。

根据辛甫森求积公式 function s=simp(x,y) if mod(length(x),2)==0 error('数据点必须为奇数') end n=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); for j=1:m-1 s2=s2+y(2*j+1); s=(y(1)+y(n)+4*s1+2*s2)*h/3; 当被积函数不是解析表示时, 比如离散数据表表示的函数 通常就用这个函数按辛甫森 公式计算积分。

二 高斯（Gauss)求积公式 各种近似求积公式都可以表示为 (11) 如何选择节点xi 和系数Ai ，使(11)计算的精度更高？ 我们不妨只考虑 而构造代数精度为3的形如 令f(x)= xk ，用(11)式计算 若对于k=0,1,...,m 都 有In = I ，而当k =m+1时，In ≠ I ，则称In 的代数精度为m 。梯 形公式代数精度为 1，辛甫森公式的代数精度为 3。 (12) G2=A1 f(x1)+ A2 f(x2) 的求积公式。 对于 f (x)=1, x, x 2, x 3，应该有 下面介绍的是取消对区间等分的限制，n 给定后同时确定 节点xi和系数Ai，使代数精度尽可能高的所谓高斯公式。 成立，依次将f(x)= 1, x, x 2, x 3代入，即可得到确定A1,A2 ,x1 ,x2 的方程组。 我们先考虑节点数为2 而使用(11)计算的积分近似值有代数 精度为3。

将 f (x) =1, x , x 2 , x 3 依次代入，得： n=2的高斯公式为： 提高精度可以通过增加节点数n ，(11)的代数精度可达到2n-1, 但增加了解高维线性方程组的难度，实用价值不大；另一方法 是将区间分小，在小区间上用G2。

将区间(a, b)作m等分，记h=(b-a)/m，xk=a+kh, k=0,1,...,m, 作变换 将x∈ [xk-1, xk]化为t∈[-1, 1] 其中 常用的高斯公式就是： 下面我们来编写M文件，应用高斯公式计算定积分。

function s=gaussinteg(f,a,b,m) s=0; h=(b-a)/m;x=a:h:b; for j=1:m 根据 定义任何一个被积函数 定义高斯积分函数M文件： function s=gaussinteg(f,a,b,m) s=0; h=(b-a)/m;x=a:h:b; for 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); s=s+feval(f,z1(j))+feval(f,z2(j)); end s=s*h/2; function y=jifenfun(x) y=exp(x).*sin(x); return 在命令窗口运行 s=gaussinteg('jifenfun', ... 0,2*pi,1000) s = -267.2458 Gauss-lobatto是改进的高斯积分方法,采取自适应求积方法

三 用MATLAB实现定积分计算： ⑴ 矩形公式与梯形公式 z1 = 0.9602 z2 = 1.0388 z11 = z12 = z3 = 0.9995 h=pi/40; x=0:h:pi/2; y=sin(x); z1=sum(y(1:20))*h, z2=sum(y(2:21))*h, z=cumsum(y)*h, z11=z(20)*h, z12=(z(21)-z(1))*h, z3=trapz(y)*h, z= 0.0062 0.3802 0.0184 0.4438 0.0368 0.5107 0.0611 0.5807 0.0911 0.6533 0.1268 0.7280 0.1678 0.8043 0.2140 0.8819 0.2650 0.9602 0.3205 1.0388

⑵ 用辛甫森公式求定积分 z4=quad('sin',0,pi/2) z5=quadl('sin',0,pi/2) z4 = 1.0000 z5 = ⑶ 高斯积分公式求定积分 前面已经保存了文件gaussinteg.m在搜索路径中，即可 直接调用。这里被积函数为内部函数，无需另外定义。 s=gaussinteg(‘sin', 0, pi/2,1000) s = 1.0000

§2 数值积分应用问题举例 a=7782.5 km b=7721.5 km 一 求卫星轨道长度 一 求卫星轨道长度 我国第一颗人造地球卫星轨道为平面椭圆曲线，近地点距 地球表面439km，远地点距地面2384km，地球半径为6371km， 求卫星轨道的长度。 椭圆的参数方程为： a=7782.5 km b=7721.5 km 椭圆长度由以下椭圆积分表示：

我们用梯形公式和辛甫生公式来求卫星轨道的长度。 t=0:pi/10:pi/2; y1=weixing(t); format long e l1=4*trapz(t,y1) l2=4*quad(‘weixing',0,pi/2,1e-6) 先创建一个函数M 文件：weixing.m function y=weixing(t) a=7782.5;b=7721.5; y=sqrt(a^2*sin(t).^2+... b^2*cos(t).^2); 运行结果： I1 = 4.870744099902405e+004 I2 = 4.870744099902417e+004 再建立一个脚本文件 tracklong.m

广义积分、二维数值积分 1、无穷区间 截断误差： blquad triplequad

上机作业 p65 3 , 4 2) ,10 请将simp.m, gaussinteg.m, 编辑好保存在自己的文件夹中，以便随时调用。 3、要求用一个脚本文件把三种公式计算的结果和计算出的精确结果用矩阵表示，在输出中能直观比较。 4、用梯形公式计算要求用不同的步长计算作比较，用辛甫森公式计算要求用不同的精度计算 作比较。设计好M文件保存输出结果便于比较。 请将simp.m, gaussinteg.m, 编辑好保存在自己的文件夹中，以便随时调用。

四 蒙特卡罗方法(用随机模拟方法计算数值积分) 1、随机投点法： 四 蒙特卡罗方法(用随机模拟方法计算数值积分) 1、随机投点法： 例、向图中边长为 1 的正方形 里随机投 n 块石头，落在四分 之一圆内的石头为 k 块，根据几何概率，则四分之一圆的面积的近似值就是k / n，即 应用：先由计算机随机产生n个点的坐标（xi , yi )(i=1,2,…n),其中xi ,yi 分别为[a，b]和[c，d]区间上的均匀分布随机数，然后记录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; end s=k/n*(b-a)*(d-c) b d a

2 均值估计法： 若随机变量X 的概率分布密度是 p(x), a x b，则 特别当X为[a, b]区间均匀分布的随机变量时，p(x)=1/(b-a), (axb) n=100000; x=rand(1,n)* pi/2; y=sin(x); z=sum(y)*pi/2/n z = 1.0010 注1、xk(k=1,2,…,n)为[a,b] 的随机数 注2、无须作随机数 yi ≤ f (xi ) 比较。 这两种用随机模拟的方式求积分近似值的方法 蒙特卡罗方法

3、蒙特卡罗方法的通用函数与调用格式 随机投点法 均值估计法 (设0≤ f(x) ≤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)<=feval(f,a+(b-a).*x(i)) k=k+1; end z=k/n*(b-a); 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); end z=z*(b-a)/n; >>z=motc1('sin',0,pi/2,1e4) >>z=motc2('sin',0,pi/2,10000) z = 1.00130915270174 z = 0.98897336735007 通常还要先设法使被积函数在[0, 1] 范围取值再应用随机投点法。

4、运用蒙特卡罗方法计算重积分 a、随机投点法 例、求半径为 的半球的体积。 function s=motc3(f,x,y,z,n) 例、求半径为 的半球的体积。 function s=motc3(f,x,y,z,n) If nargin<5 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; end s=k/n*xh*yh*zh;

b、均值估计法 计算二重积分： 设xi, yi(i =1,…,n)是相互独立，分别为[a, b]和[c,d]的随机数，判断每个点(xi, yi) 是否落在W域内，将落在 W域内的m个点记作(xk, yk) ,k=1,…,m，则 蒙特卡罗方法: 优点:计算简单,可以计算重积分; 缺点:计算量大，精度较差，结果随机

§3 数值微分 一 数值微分公式 当函数以离散数值形式给出时，用离散方法 近似地计算函数y= f (x)在某点的导数值，这就是数值微分。 一 数值微分公式 当函数以离散数值形式给出时，用离散方法 近似地计算函数y= f (x)在某点的导数值，这就是数值微分。 导数有以下三个近似公式： T y 以上三个公式分别 表示了线段AB、CA 和CB的斜率。前、 后差公式误差为O(h) 中点公式误差O(h2) y=f(x) A B 它们依次称为前差公式、后差公式和中点公式。这些公式都 明确的几何意义： C o a-h a a+h x

将区间[a, b] n等分，步长 当函数 y= f (x)在分点 上用离散数值表示为( xk , yk)，a=x0< x1 <…< xn = b 时, 函数在分点的导数值由中点公式得到三点公式： 后两个公式是由二次插值得到，目的是保持在端点处的精度O(h2) 高阶导数的近似公式一般要用插值多项式得到，下面是二阶 导数的近似公式： diff(x): [x2-x1 , x3-x2,… xn-xn-1]

已知20世纪美国人口统计数据如下，试计算表2 中这些年份的人口增长率。 二 应用实例——人口增长率 已知20世纪美国人口统计数据如下，试计算表2 中这些年份的人口增长率。 表2 20世纪美国人口统计数据 年份 1900 1910 1920 1930 1940 1950 1960 1970 1980 1990 人口(.106) 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 又已知某地区20世纪70年代的人口增长率如表3，且1970 年人口为210（百万），试估计1980年的人口。 表3 某地区20世纪70年代人口增长率数据 年份 1970 1972 1974 1976 1978 1980 年增长率(%) 0.87 0.85 0.89 0.91 0.95 1.10

1 记时刻t 的人口为x(t)，人口相对增长率为 ，记 1900—1990年的人口依次为 xk ，(k=0,1,...,n)，年增长率为rk 。 由三点公式可以得到： r = 0.0220 0.0166 0.0146 0.0102 0.0104 0.0158 0.0149 0.0116 0.0105 x=[76.0 92.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); end r(10)=(x(8)-4*x(9)+3*x(10))/20/x(10); r

2 某地区20世纪70年代人口增长率如下表： 人口增长满足微分方程 和初始条件为x(0)= x0。 其解为 由于增长率r(t)为离散数据，故 2 某地区20世纪70年代人口增长率如下表： 表3 某地区20世纪70年代人口增长率数据 年份 1970 1972 1974 1976 1978 1980 年增长率(%) 0.87 0.85 0.89 0.91 0.95 1.10 人口增长满足微分方程 和初始条件为x(0)= x0。 其解为 由于增长率r(t)为离散数据，故 用数值积分计算： x0=210; r=[0.87 0.85 0.89 0.91 0.95 1.10]/100; y=trapz(r)*2; %步长为2 x=x0*exp(y) x = 230.1676

实验内容与要求 p93 5.2 ２) b，d 5) 10) 10) 概率密度函数是 其中σx, σy分别为x方向和y方向的均方差,ρ为偏差在x和y之间的相关系数以100m为一个单位，目标区x^2+y^2<=1 采用蒙特卡罗方法之均值估计法，参考P91的4.2射击命中概率 11)用蒙特卡罗法计算重积分。(思考)