ch4_数值计算2.doc

PAGE  PAGE 32 数值计算 本章主要讲的函数有 1、 近似数值极限及导数的函数有 diff gradient 2 、 数值求和与近似数值积分函数有 sum sumsum trapz cumtranpz 3、 计算精度可控的数值积分 integral 4、 函数极值的数值求解 fminbnd fminserch 5、 常微分方程的数值求解 ode45 6、 矩阵和代数方程 求矩阵的秩：rank 迹：trace 行列式：det 阶梯矩阵变换：rref A矩阵零空间的全部正交基：null A矩阵值空间的全部正交基：orth A矩阵的特征值，特征向量分解：eig 7、 一般代数方程的解 一元函数零点指令：fzero 解非线性方程组：fsolve 数值微积分 （1）重视有限精度浮点表示的离散本质，不要贸然自行编写数值计算程序进行求导和求极限运算 （2）数值导数一定要求取，自变量的增量大于原数据精度的10倍以上。 （3）解极限值，积分，微分方程数值计算时，要尽量使用MATLAB提供的指令，严格遵守指令的使用规则。 特别说明： MATLAB没有提供数值求导和极限的指令。 其中概念法求差分 dx=diff(X)求差分 FX=gradient(F) [FX,FY]=gradient(F) 近似数值极限及导数 在数值法近似求极限和理论有时可能不一致（举例说明）,数值法求极限轻易不用 【例4.1-1】设试用机器零阈值eps替代理论0计算极限 理论分析： %不可信的极限的数值近似计算 x=eps; L1=(1-cos(2*x))/(x*sin(x)), % L2=sin(x)/x, % L1 = 0 L2 = 1 %可信的极限的符号计算 syms t f1=(1-cos(2*t))/(t*sin(t)); f2=sin(t)/t; Ls1=limit(f1,t,0) Ls2=limit(f2,t,0) Ls1 = 2 Ls2 = 1 【例4.1-2】已知x=sin(t)，求该函数在区间[0 2pi]中的近似导函数。 本例考察的是，在将连续信号进行离散化时，采样间隔影响导函数。 （1）计算数值导数时，自变量的额增量取得过小（在eps数量级） d=pi/100; t=0:d:2*pi; x=sin(t); dt=5*eps; %增量为eps x_eps=sin(t+dt); dxdt_eps=(x_eps-x)/dt; %以dt=5*eps为增量算得的数值导数 plot(t,x,LineWidth,5) hold on plot(t,dxdt_eps) hold off legend(x(t),dx/dt) xlabel(t) 图 4.1-1 增量过小引起有效数字严重丢失后的毛刺曲线 （2）计算数值导数时，自变量的增量取得适当 x_d=sin(t+d); dxdt_d=(x_d-x)/d; %以d=pi/100为增量得的数值导数 plot(t,x,LineWidth,5) hold on plot(t,dxdt_d) hold off legend(x(t),dx/dt) xlabel(t) 图 4.1-2 增量适当所得导函数比较光滑 现象差距合理解释 如果dt步长太小，f(t+dt)与f(t)数值十分接近，它们的高维有效数字完全相同。f(t+dt)与f(t)差高位有效数字消失，导致精度急剧下降。 【例4.1-3】已知x=sin(t),采用diff和gradient计算该函数在区间[0 2*pi]中的近似导数 clf d=pi/100; % t=0:d:2*pi; x=sin(t); dxdt_diff=diff(x)/d; % dxdt_grad=gradient(x)/d; % subplot(1,2,1) plot(t,x,b) hold on plot(t,dxdt_grad,m,LineWidth,8) plot(t(1:end-1),dxdt_diff,.k,MarkerSize,8) axis([0,2*pi,-1.1,1.1]) title([0, 2\pi]) legend(x(t),dxdt_{grad},dxdt_{diff},Location,North) xlabel(t),box off hold off subplot(1,2,2) kk=(length(t)-10):length(t);%t hold on plot(t(kk),dxdt_grad(kk),om,MarkerSize,8) plot(t(kk-1),dxdt_diff(kk-1),.k,MarkerSiz