《非数与数值微积分》PPT课件.ppt

第七讲 非数与 数值微积分,本讲教学目标 掌握非数的概念和性质 掌握空数组的概念和创建方法 掌握数值微分、梯度和法线 掌握常见的数值积分方法 了解二重积分的方法,7.1 “非数”和“空数组” 在进行科学计算时经常遇到除零的现象，这时就会出现错误信息。 MATLAB中规定出现上述情况不出致命错误信息，仅在输出中进行警告并指出这是非数。 我们可以巧妙地利用非数对图形进行镂空、裁减等操作，从而得到很好的可视化效果。,1.1.1 “非数” MATLAB 规定，形如 都会产生非数(Not a Number)，用NaN或nan 表示。 根据 IEEE 数学规范，NaN具有以下性质： 只要NaN 参加四则运算或函数作用于非数其结果也一定是非数，即具有传递性。 非数没有大小的概念，因此不能进行逻辑运算也不能比较两个非数的大小。,例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 = NaN 1）非数的传递性 0*a,sin(a) ans = NaN ans = NaN,2）非数的不可比较性 a=nan % 非数不能比较 ans = 0 a=nan % 非数不能进行逻辑运算 bc % 两个非数不能比较大小 ans = 1 ans = 0 3）非数的属性判断 class(a) % 数据类型归属 isnan(a) % 唯一正确判断非数的指令 ans = double ans = 1,从结果可知非数也是一个数，数据类型是双精度的。,7.1.2 “空”数组 当数组的某些维数为0时，称该数组为空数组。 在复杂的计算中有时会产生空数组，也可以利用空数组对其他数组进行赋值。 空数组与全零元素数组不是一回事 空数组确实是数组，用which、who、whos可以看出。 可以用命令 isempty 来判断一个数组是否为空数组。,例2 创建空数组。 A = rand(2,2,2); % 直接赋值法 A(:,:,:) = ; ndims(A) size(A) ans = 3 ans = 0 2 2 可看到 A 仍然是 3 维数组，但其中有一个维数是 0，则是空数组。,(2) 用函数创建空数组 A = ones(25,0,4); if A S1=是数组 else S0=是空数组 end S0 = 是空数组,例3 “空”数组的属性。 class(A) % 显示数据类别 isnumeric(A) % 是数值数组类吗 isempty(A) % 是否为空数组？ size(A) ans =double ans = 1 ans = 1 ans = 25 0 4,从结果看，数组A 是一个双精度的空数组，具有维数2504。,7.2 数值微积分 7.2.1 数值微分 在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函数diff，其格式为： DX=diff(X)：计算向量X的向前差分 DX(i)=X(i+1)-X(i)，i=1,2,n-1。 DX=diff(X,n)：计算X的n阶向前差分 diff(X,2)=diff(diff(X)。 DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(缺省)，按列计算；dim=2，按行计算。,例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(发展加速度),差分方法应用十分广泛，如在经济数据中一阶差分表示发展速度，二阶差分表示经济发展的加速度。,7.2.2 梯度 在MATLAB中，使用gradient函数求解近似梯度，其格式为： fx, fy = gradient(f) 返回矩阵 f 的数值梯度，fx相当于df/dx，即在x方向（列）的差分值。 fx, fy = gradient(f, h) 求二元函数的梯度，使用数量 h 作为各个方向的间隔点。 fx, fy, fz, = gradient(f, h1, h2, ) 做相似工作。,例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),7.2.3 法线 计算与绘制法线的命令为： NX, NY, NZ = surfnorm( X, Y, Z ) 方向法线绘图命令的基本格式为： quiver3(Z,U,V,W) quiver3(X,Y,Z,U,V,W) quiver3(.,scale) 其中：X,Y,Z 为曲线上的点；U,V,W 为每一点上切面的法线方向；scale 为实数值，表示方向箭头的指向。,例6 求函数 的曲面图，且作出每点的切面法向量图形。 X,Y = meshgrid(-2:0.25:2,-1:0.2:1); Z = X.* exp(-X.2 - Y.2); U,V,W = surfnorm(X,Y,Z); % 计算法向量分量 quiver3(X,Y,Z,U,V,W,0.5); % 作法向量的图形 hold on surf(X,Y,Z);colormap hsv view(-35,45) axis (-2 2 -1 1 -.6 .6) hold off,7.2.4 数值积分 求解定积分的数值方法多种多样，如下都是经常采用的方法。 简单的梯形法 辛普生(Simpson)法 牛顿柯特斯(Newton-Cotes)法 其基本思想都是将积分区间a,b分成n个子区间xi,xi+1，i=1,2,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。,1. 变步长辛普生法 调用格式为： q = quad(fun,a,b) q = quad(fun,a,b,tol) q,n=quad(fun,a,b,tol,trace) 其中fun为被积函数；a, b是积分上下限，tol为给定的误差限，缺省取0.001；trace若取0则不展现积分过程，缺省取0；取非0则展现。返回参数q即定积分值，n为被积函数的调用次数。,例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,例8 对空间曲线求积分。 设空间曲线的参数方程为： x = t.*sin(t), y = t.*cos(t), z=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),计算结果为： len = 8.5191,2. 牛顿柯特斯法 调用格式为： I,n=quad8(fun,a,b,tol,trace) 其中参数的含义和quad函数相似，只是tol的缺省值取10-6。 该函数可以更精确地求出定积分的值，且一般情况下函数调用的步数明显小于quad函数，从而保证能以更高的效率求出所需的定积分值。,例9 分别用函数quad和quad8求定积分的近似值，并在相同的精度下比较函数调用次数。 format long; %调用函数quad求定积分 fx=inline(exp(-x); I,n=quad(fx,1,2.5,1e-10) I = 0.28579444254766 n = 65, format long; %调用函数quad8 fx=inline(exp(-x); I,n=quad8(fx,1,2.5,1e-10) I = 0.28579444254754 n = 33,3. 被积函数由一个表格定义 用trapz(X,Y)函数求解，其中向量X,Y定义函数关系Y=f(X)。 例10 用trapz函数计算定积分 。 X=1:0.01:2.5; Y=exp(-X); % 生成函数关系数据向量 trapz(X,Y) ans = 0.28579682416393,4. 二重积分 与一重积分的原理一样，其语法为：