第7讲：非数与数值微积分

1、第七讲第七讲 非数与非数与 数值微积分数值微积分2n本讲教学目标本讲教学目标 掌握非数的概念和性质掌握非数的概念和性质 掌握空数组的概念和创建方法掌握空数组的概念和创建方法 掌握数值微分、梯度和法线掌握数值微分、梯度和法线 掌握常见的数值积分方法掌握常见的数值积分方法 了解二重积分的方法了解二重积分的方法3n7.1 “非数非数”和和“空数组空数组”在进行科学计算时经常遇到除零的现象，这时在进行科学计算时经常遇到除零的现象，这时就会出现错误信息。就会出现错误信息。matlab中规定出现上述情况不出致命错误信中规定出现上述情况不出致命错误信息，仅在输出中进行警告并指出这是非数。息，仅在输出中进行警

2、告并指出这是非数。我们可以巧妙地利用非数对图形进行镂空、裁我们可以巧妙地利用非数对图形进行镂空、裁减等操作，从而得到很好的可视化效果。减等操作，从而得到很好的可视化效果。4n1.1.1 “非数非数”matlab 规定，形如规定，形如 都会产生非数都会产生非数(not a number)，用，用nan或或nan 表示。表示。根据根据 ieee 数学规范，数学规范，nan具有以下性质：具有以下性质：n只要只要nan 参加四则运算或函数作用于非数其结果也参加四则运算或函数作用于非数其结果也一定是非数，即具有传递性。一定是非数，即具有传递性。n非数没有大小的概念，因此不能进行逻辑运算也不非数没有大小的

3、概念，因此不能进行逻辑运算也不能比较两个非数的大小。能比较两个非数的大小。0,005n例例1 非数的产生和性质演示。非数的产生和性质演示。 a=0/0,b=0*log(0),c=inf/inf ,d=sqrt(nan) warning: divide by zero. in d:zjjpaperch23.m at line 62 a = nan; b = nan; c = nan; d = nan1）非数的传递性）非数的传递性 0*a,sin(a) ans = nan ans = nan6n2）非数的不可比较性）非数的不可比较性 a=nan % 非数不能比较非数不能比较 ans = 0 a=n

4、an % 非数不能进行逻辑运算非数不能进行逻辑运算 bc % 两个非数不能比较大小两个非数不能比较大小 ans = 1 ans = 0n3）非数的属性判断）非数的属性判断 class(a)% 数据类型归属数据类型归属 isnan(a)% 唯一正确判断非数的指令唯一正确判断非数的指令 ans = double ans = 1 n从结果可知非数也从结果可知非数也是一个数，数据类是一个数，数据类型是双精度的。型是双精度的。7n7.1.2 “空空”数组数组当数组的某些维数为当数组的某些维数为0时，称该数组为空数组时，称该数组为空数组。在复杂的计算中有时会产生空数组，也可以利在复杂的计算中有时会产生空数

5、组，也可以利用空数组对其他数组进行赋值。用空数组对其他数组进行赋值。空数组与全零元素数组不是一回事空数组与全零元素数组不是一回事空数组确实是数组，用空数组确实是数组，用which、who、whos可可以看出。以看出。可以用命令可以用命令 isempty 来判断一个数组是否为空来判断一个数组是否为空数组。数组。8n例例2 创建空数组。创建空数组。 a = rand(2,2,2); % 直接赋值法直接赋值法 a(:,:,:) = ; ndims(a) size(a) ans = 3 ans = 0 2 2可看到可看到 a 仍然是仍然是 3 维数组，但其中有一个维数维数组，但其中有一个维数是是 0，

6、则是空数组。，则是空数组。n(2) 用函数创建空数组用函数创建空数组 a = ones(25,0,4); if a s1=是数组是数组 else s0=是空数组是空数组 ends0 = 是空数组是空数组9n例例3 “空空”数组的属性。数组的属性。 class(a) % 显示数据类别显示数据类别 isnumeric(a) % 是数值数组类吗是数值数组类吗 isempty(a) % 是否为空数组？是否为空数组？ size(a) ans =double ans = 1 ans = 1 ans = 25 0 4n从结果看，数组从结果看，数组a 是一是一个双精度的空数组，具个双精度的空数组，具有维数有维

7、数2504。10n7.2 数值微积分数值微积分7.2.1 数值微分数值微分n在在matlab中，没有直接提供求数值导数的函数，中，没有直接提供求数值导数的函数，只有计算向前差分的函数只有计算向前差分的函数diff，其格式为：，其格式为：ndx=diff(x)：计算向量计算向量x的向前差分的向前差分 dx(i)=x(i+1)-x(i)，i=1,2,n-1。ndx=diff(x,n)：计算计算x的的n阶向前差分阶向前差分 diff(x,2)=diff(diff(x)。ndx=diff(a,n,dim)：计算矩阵计算矩阵a的的n阶差分，阶差分，dim=1时时(缺省缺省)，按列计算；，按列计算；dim

8、=2，按行计算。，按行计算。11n例例4 对我国对我国5298年国民经济数据年国民经济数据gdp发展趋势发展趋势进行分析，求其一、二阶差分，并作图。进行分析，求其一、二阶差分，并作图。 a=load(e:datagdp.txt) b=a(:,2) b1=diff(b); b2=diff(b,2) subplot(1,3,1),plot(b),title(发展趋势发展趋势) subplot(1,3,2),plot(b1),title(发展速度发展速度) subplot(1,3,3),plot(b2),title(发展加速度发展加速度)n差分方法应用十分广差分方法应用十分广泛，如在经济数据中泛，如

9、在经济数据中一阶差分表示发展速一阶差分表示发展速度，二阶差分表示经度，二阶差分表示经济发展的加速度。济发展的加速度。12n7.2.2 梯度梯度在在matlab中，使用中，使用gradient函数求解近似梯函数求解近似梯度，其格式为：度，其格式为：nfx, fy = gradient(f) 返回矩阵返回矩阵 f 的数值梯度，的数值梯度，fx相当于相当于df/dx，即在，即在x方向（列）的差分值。方向（列）的差分值。nfx, fy = gradient(f, h) 求二元函数的梯度，求二元函数的梯度，使用数量使用数量 h 作为各个方向的间隔点。作为各个方向的间隔点。nfx, fy, fz, = g

10、radient(f, h1, h2, ) 做相做相似工作。似工作。13n例例5 对函数对函数 作梯度图。作梯度图。 v = -2:0.2:2; x,y = meshgrid(v); z = x .* exp(-(x.2 + y.2); subplot(1,2,1);surf(x,y,z), subplot(1,2,2); px,py = gradient(z,.2,.2); contour(x,y,z), hold on, quiver(x,y,px,py)22xyxe14n7.2.3 法线法线计算与绘制法线的命令为：计算与绘制法线的命令为：n nx, ny, nz = surfnorm( x

11、, y, z )方向法线绘图命令的基本格式为：方向法线绘图命令的基本格式为：nquiver3(z,u,v,w)nquiver3(x,y,z,u,v,w)nquiver3(.,scale)n其中：其中：x,y,z 为曲线上的点；为曲线上的点；u,v,w 为每一点上切为每一点上切面的法线方向；面的法线方向；scale 为实数值，表示方向箭头的为实数值，表示方向箭头的指向。指向。15n例例6 求函数求函数 的曲面图，且作出每点的切的曲面图，且作出每点的切面法向量图形。面法向量图形。x,y = meshgrid(-2:0.25:2,-1:0.2:1);z = x.* exp(-x.2 - y.2);u

12、,v,w = surfnorm(x,y,z); % 计算法向量分量计算法向量分量quiver3(x,y,z,u,v,w,0.5); % 作法向量的图形作法向量的图形hold onsurf(x,y,z);colormap hsvview(-35,45)axis (-2 2 -1 1 -.6 .6)hold off22xyxe16n7.2.4 数值积分数值积分求解定积分的数值方法多种多样，如下都是经求解定积分的数值方法多种多样，如下都是经常采用的方法。常采用的方法。n简单的梯形法简单的梯形法n辛普生辛普生(simpson) 法法n牛顿柯特斯牛顿柯特斯(newton-cotes)法法其基本思想都是将

13、积分区间其基本思想都是将积分区间a,b分成分成n个子区个子区间间xi,xi+1，i=1,2,n，其中，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。这样求定积分问题就分解为求和问题。17n1. 变步长辛普生法变步长辛普生法调用格式为：调用格式为： q = quad(fun,a,b) q = quad(fun,a,b,tol) q,n=quad(fun,a,b,tol,trace)n其中其中fun为被积函数；为被积函数；a, b是积分上下限，是积分上下限，tol为为给定的误差限，缺省取给定的误差限，缺省取0.001；trace若取若取0则不则不展现积分过程，缺省取展现积分过程，缺

14、省取0；取非；取非0则展现。返回则展现。返回参数参数q即定积分值，即定积分值，n为被积函数的调用次数。为被积函数的调用次数。18n例例7 计算积分计算积分 q = quad(1./(1+25*x.2),-1,1) ans = 0.5494% 也可以使用内联函数定义被积函数也可以使用内联函数定义被积函数 f = inline(1./(1+25*x.2); q = quad(f,-1,1); ans = 0.5494 1122511dxx19n例例8 对空间曲线求积分。对空间曲线求积分。设空间曲线的参数方程为：设空间曲线的参数方程为： x = t.*sin(t), y = t.*cos(t), z

15、=t t = 0:pi/50:10*pi; plot3(t.*sin(t),t.*cos(t),t,linewidth,2) f = inline(sqrt(4*(t.*sin(t).2 + (t.*cos(t).2 + 1); len = quad(f,0,pi)n计算结果为：计算结果为： len = 8.519120n2. 牛顿柯特斯法牛顿柯特斯法调用格式为：调用格式为： i,n=quad8(fun,a,b,tol,trace)n其中参数的含义和其中参数的含义和quad函数相似，只是函数相似，只是tol的缺省的缺省值取值取10-6。 该函数可以更精确地求出定积分的值，且一般该函数可以更精确

16、地求出定积分的值，且一般情况下函数调用的步数明显小于情况下函数调用的步数明显小于quad函数，从函数，从而保证能以更高的效率求出所需的定积分值。而保证能以更高的效率求出所需的定积分值。21n例例9 分别用函数分别用函数quad和和quad8求定积分的近似值，求定积分的近似值，并在相同的精度下比较函数调用次数。并在相同的精度下比较函数调用次数。 format long; %调用函数调用函数quad求定积分求定积分 fx=inline(exp(-x); i,n=quad(fx,1,2.5,1e-10) i = 0.28579444254766 n = 65n format long; %调用函数调用函数quad8n fx=inline(exp(-x);n i,n=quad8(fx,1,2.5,1e-10) i = 0.28579444254754 n = 3322n3. 被积函数由一个表格定义被积函数由一个表格定义用用trapz(x,y)函数求解，其中向量函数求解，其中向量x,y定义函数定义函数关系关系y=f(x)。n例例10 用用trapz函数计算定积分函数计算定积分 。 x=1:0.01:2.5; y=exp(-x); % 生成函数关系数据向量生成函数关系数据向量 trapz(x,y) ans = 0.285796824163935 . 21dxex23n4. 二重积分二重积分