[理学]第8讲 数值积分与微分.ppt

人口增长率 已知20世纪美国人口统计数据如下，试计算表中这些年份的人口增长率 解：若记时刻t的人口为x(t)，则人口(相对)增长率为r(t)=(dx/dt)/x(t)，表示每年人口增长的比例，对于题目给出的人口数据，记1990年为k=0, 1910,1920,…,1990年依次为k=1,2,…,9，相应地，人口记为xk，年增长率为rk，用数值微分的三点公式，即 数值积分 定积分的定义： 2.1 梯形公式和辛普森公式 2.1.1 公式的导出 曲边梯形的面积 若将二者平均，则每个小区间上的小矩形变为小梯形（图中粗实线），整个区间上的结果为 求m段之和就得到整个区间上的近似积分公式 2.1.2 误差估计和收敛性 梯形公式在小区间[xk,xk+1]?上是用线性插值函数T(x)代替f(x)，容易得到 进一步，记 类似地，可以求出辛普森求积公式（5）的误差 2.1.3 步长的自动选择 2.2 高斯求积公式 2.2.1 代数精度 从（2）-（5）式我们看到，不论哪个近似求积公式都可以写成如下形式： 代数精度是衡量求积公式精确程度的一种指标，它用密函数作为被积函数f，以近似积分与精确值是否相等作为精度的度量标准，有如下定义： 设 2.2.2 高斯公式 区间(a,b)经过适当的变量代换可以化为(-1,1)，因而只需计算 为了提高精度，可采用两种方法：一是增加n；但要解形如（13）的、更复杂的非线性方程组，n较大时实用价值不大；另一个是将区间(a,b)分小，在小区间上用G2，具体做法如下：将区间（a，b）m等分，记 用G2近似Ik，就有 这就是常用的高斯公式。 2.3 蒙特卡罗法 举例：在一个1/4单位圆中，如果向图中边长为1的正方形内随机投n块小石头，当n很大时小石头会均匀分布在正方形中，数一下落在1/4圆内的石头，假定有k个，那么k/n就看作1/4单位圆面积 的近似值，于是有 ，显然，这可以看作近似计算 的一种方法。 2.3.1 随机投点法 从概率论的观点看，记投点的坐标为 于是当 这里n是（二维）（0,1）随机数（xi,yi）的总数，k是其中满足yif(xi)的数目，(0,1)随机数用计算机可以方便地产生。 2.3.2 均值估计法 若随机变量X的概率分布密度p(x),a≤x≤b,则 2.4 用Matlab作数值积分 sum(x)：输入数组x，输出为x的和，可用于按矩形公式（2）、（3）计算积分。 cumsum(x)：输入数组x，输出数组为x的依次累加和，也可用于按矩形公式（2）、（3）计算积分。 trapz(x)：输入数组x，输出为按梯形公式（4）计算的x的积分（单位步长）。 Trapz(x,y)：输入x,y为同长度的数组，输出y对x的积分（按梯形公式计算，但步长不一定相等）。 quad(‘fun’,a,b)：用辛普森（2阶）公式计算以fun.m文件命名的函数（或库函数如‘sin’,’log’）在区间(a,b)上的积分，自动选择步长，相对误差为10-3，输出积分值。 quad(‘fun’,a,b,tol)：同上，但指定相对误差为tol。 Quad8(‘fun’,a,b,tol)：用辛普森（8阶）公式计算，在同样精度下需要节点数较小。 rand(1,n)：产生n个(0,1)随机数，用于蒙特卡罗法。 例 用Matlab计算 解：用以下几种公式计算： （1）矩形公式和梯形公式 将 10等分，步长为 (2)辛普森公式 z4=quad(‘sin’,0,pi/2) 即得z4=1.0000，这显然更方便和精确，但是quad和quad8不能像trapz那样，对用数值给出的x,y数组作积分。 （3）蒙特卡罗法 用均值估计法，编程如下： n=10000; x=rand(1,n); y=sin(x.*pi/2); z=sum(y)*pi/2/n 所得的结果是随机的。 §3 数值微分 数值积分是用离散方法近似地计算函数y=f(x)在某点x=a的导数值。根据导数定义，可用差商近似微商（导数），有 如果将二者平均一下，可得： 解：此题无法用解析法求解，用蒙特卡罗法求解。 作变换x=au,y=bv,且以100(m)为1单位，即σx= σy=1,a=1.2, b=0.8.于是 程序： a=1.2;b=0.8;m=0;z=0;n=100000; for i=1:n x=rand(1,2) y=0; if x(1)^2+x(2)^2=1 y=exp(-0.5*(a^2*x(1)^2+b^2*x(2)^2)); z=z+y; m=m+1; end end