第4MATLAB的数值计算.ppt

第4MATLAB的数值计算 Still waters run deep.流静水深流静水深,人静心深人静心深 Where there is life,there is hope。有生命必有希望。有生命必有希望4.1数值微积分数值微积分4.1.1 近似数值极限及导数近似数值极限及导数4.1.2 数值求和与近似数值积分数值求和与近似数值积分4.1.3 计算精度可控的数值积分计算精度可控的数值积分4.1.4 函数极值的数值求解函数极值的数值求解4.1.5 常微分方程的数值解常微分方程的数值解12/6/2022第2页u在在MATLAB数值计算中，既没有专门的求极限指令，数值计算中，既没有专门的求极限指令，也没有专门的求导指令。但也没有专门的求导指令。但MATLAB提供了与提供了与“求导求导”概念有关的概念有关的“求差分求差分”指令。指令。udx=diff(X)%计算向量计算向量X的前向差分的前向差分uFX=gradient(F)%求一元求一元(函数函数)梯度梯度uFX,FY=gradient(F)%求二元（函数）梯度求二元（函数）梯度u对对diff而言，当而言，当X是向量时，是向量时，dx=X(2:n)-X(1:n-1);当当X是矩阵时，是矩阵时，dx=X(2:n,:)-X(1:n-1,:)。dx的长度比的长度比x的长度少的长度少1个元素个元素。4.1.1 近似数值极限及导数近似数值极限及导数12/6/2022第3页udx=diff(X)%计算向量计算向量X的前向差分的前向差分uFX=gradient(F)%求一元求一元(函数函数)梯度梯度uFX,FY=gradient(F)%求二元（函数）梯度求二元（函数）梯度u对对gradient而言，当而言，当F是向量时，是向量时，FX(1)=F(2)-F(1),FX(2:end-1)=(F(3:end)-F(1:end-2)/2,FX(end)=F(end)-F(end-1);FX长度与长度与F的长度相同的长度相同u当当F是矩阵时，是矩阵时，FX,FY是与是与F同样大小的矩阵。同样大小的矩阵。FX的的每行对应每行对应F相应行元素间的梯度相应行元素间的梯度;FY的每列对应的每列对应F相应相应列元素间的梯度列元素间的梯度;4.1.1 近似数值极限及导数近似数值极限及导数12/6/2022第4页数值极限和导数的应用应十分谨慎数值极限和导数的应用应十分谨慎x=eps;L1=(1-cos(2*x)/(x*sin(x),L2=sin(x)/x,L1=0L2=1syms tf1=(1-cos(2*t)/(t*sin(t);f2=sin(t)/t;Ls1=limit(f1,t,0)Ls2=limit(f2,t,0)Ls1=2Ls2=1 x=pi/1000;%可得到正确结果可得到正确结果12/6/2022第5页数值极限和导数的应用应十分谨慎数值极限和导数的应用应十分谨慎%（1）自变量的增量取得过小（）自变量的增量取得过小（eps数量级）数量级）d=pi/100;t=0:d:2*pi;x=sin(t);dt=5*eps;x_eps=sin(t+dt);dxdt_eps=(x_eps-x)/dt;plot(t,x,LineWidth,5)hold onplot(t,dxdt_eps)hold offlegend(x(t),dx/dt)xlabel(t)【例【例4.1-2】已知】已知,求求该该函数在区函数在区间间 中的近似中的近似导导函数。函数。数数值导值导数受数受计计算中有限精度困算中有限精度困扰扰，当增量，当增量dt过过小小时时，f(t+dt)与与f(t)的数的数值值十分接近，高位有效数字完全相十分接近，高位有效数字完全相同，同，df=f(t+dt)-f(t)造成高位有效数字消失，精度造成高位有效数字消失，精度急急剧变剧变差。差。12/6/2022第6页数值极限和导数的应用应十分谨慎数值极限和导数的应用应十分谨慎%（2）自变量的增量取得适当）自变量的增量取得适当x_d=sin(t+d);dxdt_d=(x_d-x)/d;plot(t,x,LineWidth,5)hold onplot(t,dxdt_d)hold offlegend(x(t),dx/dt)xlabel(t)【例【例4.1-2】已知】已知,求求该该函数在区函数在区间间 中的近似中的近似导导函数。函数。d=pi/100;12/6/2022第7页d=pi/100;t=0:d:2*pi;x=sin(t);dxdt_diff=diff(x)/d;dxdt_grad=gradient(x)/d;【例【例4.1-3】已知】已知 采用采用diff和和gradient计算该计算该函数在区间函数在区间 中的近似导函数。中的近似导函数。subplot(1,2,1);plot(t,x,b);hold onplot(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,2pi)legend(x(t),dxdt_grad,dxdt_diff,Location,North)xlabel(t),box off;hold offsubplot(1,2,2)kk=(length(t)-10):length(t);hold on;plot(t(kk),dxdt_grad(kk),om,MarkerSize,8)plot(t(kk-1),dxdt_diff(kk-1),.k,MarkerSize,8)title(end-10,end)legend(dxdt_grad,dxdt_diff,Location,SouthEast)xlabel(t),box off;hold off 宏观上宏观上,diff和和gradient结果大致相同结果大致相同细节上细节上,diff和和gradient数值有差异，数值有差异，diff没有给出最后一点导数没有给出最后一点导数12/6/2022第8页Sx=sum(X)%sum按列向求和得按列向求和得(1n)数组数组Scs=cumsum(X)%沿沿X列向求列向求累计和累计和,仍是数组仍是数组,第第(i,k)个元素是个元素是X数组第数组第k列前列前i个元素的和。最后一行等于个元素的和。最后一行等于Sx4.1.2 数值求和与近似数值积分数值求和与近似数值积分St=trapz(t,X)或或 St=dt*trapz(X)%梯形法求積分梯形法求積分Sct=cumtrapz(t,X)或或 Sct=dt*cumtrapz(X)%梯形法沿列向求梯形法沿列向求X关于关于x的的累计积分累计积分,最后一个值等于最后一个值等于StS=dt*sum(X),S=sum(t,X)%近似矩形法求积分近似矩形法求积分Scs=cumsum(t,X)=dt*cumsum(X)12/6/2022第9页clear;d=pi/8;t=0:d:pi/2;y=0.2+sin(t);s=sum(y);s_sa=d*s;%s_sa=sum(t,y),近似矩形法积分近似矩形法积分s_ta=trapz(t,y);%梯形法积分梯形法积分【例【例4.1-4】求积分】求积分,其中其中disp(sum求得积分求得积分,blanks(3),trapz求得积分求得积分)disp(s_sa,s_ta)t2=t,t(end)+d;y2=y,nan;stairs(t2,y2,:k);hold onplot(t,y,r,LineWidth,3)h=stem(t,y,LineWidth,2);set(h(1),MarkerSize,10)axis(0,pi/2+d,0,1.5)hold off;shg sum求得积分求得积分 trapz求得积分求得积分 1.5762 1.301312/6/2022第10页4.1.3 计算精度可控的数值积分计算精度可控的数值积分一重积分一重积分(quadrature精度可控精度可控)：S1=quad(fun,a,b,tol)%自适应自适应Simpson法法S2=quadl(fun,a,b,tol)%自适应罗巴托自适应罗巴托 Lobatlo法法二重积分二重积分(精度可控精度可控)：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重积分三重积分(精度可控精度可控)：S4=triplequad(fun,xmin,xmax,ymin,ymax,zmix,zmax,tol)fun:可为字符串可为字符串,内联对象内联对象,匿名函数匿名函数,函数句柄函数句柄,注意乘注意乘除法运算一定加点除法运算一定加点.用数组运算用数组运算 tol:默认积分的绝对精度为默认积分的绝对精度为10-612/6/2022第11页（1）syms xIsym=vpa(int(exp(-x2),x,0,1)%x.2数组乘方亦可数组乘方亦可Isym=0.74682413281242702539946743613185 例例4.1-5：求：求 .%（2）梯形法积分梯形法积分format longd=0.001;x=0:d:1;Itrapz=d*trapz(exp(-x.*x)%x.*x必须为数组乘必须为数组乘Itrapz=0.746824071499185 （3）fx=exp(-x.2);%一定用数组乘一定用数组乘Ic=quad(fx,0,1,1e-8)Ic=0.746824132854452%实际精度控制到实际精度控制到1e-1012/6/2022第12页（1）符号计算法）符号计算法syms x ys=vpa(int(int(xy,x,0,1),y,1,2)Warning:Explicit integral could not be found.s=0.40546510810816438197801311546435 例例4.1-6：求：求 .（2）数值积分法）数值积分法format longs_n=dblquad(x,y)x.y,0,1,1,2)%匿名函数表示被积函数匿名函数表示被积函数s_n=0.405466267243508 s_n=dblquad(x.y,0,1,1,2)%字符串表示被积函数字符串表示被积函数s_n=dblquad(inline(x.y),0,1,1,2)%内联函数表示被积函数内联函数表示被积函数%一定为数组乘一定为数组乘12/6/2022第13页4.1.4 函数极函数极值值的数的数值值求解求解 x,fval,exitflag,output=fminbnd(fun,x1,x2,options)%一元函数一元函数在区间在区间bound(x1,x2)中极小值中极小值 x,fval:极值点坐标和对应目标函数极值极值点坐标和对应目标函数极值x,fval,exitflag,output=fminsearch(fun,x0,options)%单纯形法求搜索起点单纯形法求搜索起点x0附近附近多元函数多元函数极值点极值点%x每列代表一个候选极值点，各列按每列代表一个候选极值点，各列按目标函数极小值递增顺序目标函数极小值递增顺序%x(:,1)对应的目标函数极小值点由对应的目标函数极小值点由fval给出给出 fun:字符串字符串,内联函数内联函数,匿名函数匿名函数,函数句柄函数句柄,注意乘除法注意乘除法运算一定加点运算一定加点.用数组运算用数组运算 options:配置优化参数配置优化参数,可略可略 exitflag:给出大于给出大于0的数的数,则成功搜索到极值点则成功搜索到极值点 output:给出具体的优化算法和迭代次数给出具体的优化算法和迭代次数12/6/2022第14页【例【例4.1-7】已知】已知 ，在在-10 x1010区区间间，求函数的最小求函数的最小值值。（1）用）用“导数为零法导数为零法”求极值点求极值点syms xy=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);yd=diff(y,x);%求导函数求导函数xs0=solve(yd,x)%求导函数为零的根，即极值点求导函数为零的根，即极值点yd_xs0=vpa(subs(yd,x,xs0)%验证极值点处导函数为零验证极值点处导函数为零y_xs0=vpa(subs(y,x,xs0)%求极值点处极值求极值点处极值xs0=matrix(0.050838341410271656880659496266968)yd_xs0=2.2958874039497802890014385492622*10(-41)y_xs0=-0.001263317776974196724544154344118 无法判断是否最小值无法判断是否最小值12/6/2022第15页【例【例4.1-7】已知】已知 ，在在-10 x1010区区间间，求函数的最小求函数的最小值值。（2）采用优化算法求极小值）采用优化算法求极小值x1=-10;x2=10;yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);xn0,fval,exitflag,output=fminbnd(yx,x1,x2)xn0=2.514797840754235fval=%比比“导数为零法导数为零法”求得的极值更小求得的极值更小,更可能是最小值更可能是最小值 -0.499312445280039exitflag=1output=iterations:13 funcCount:14 algorithm:golden section search,parabolic interpolation message:1x112 char 12/6/2022第16页【例【例4.1-7】已知】已知 ，在在-10 x1010区区间间，求函数的最小求函数的最小值值。（4）据图形观察，重设）据图形观察，重设fminbnd的搜索区间的搜索区间x11=6;x2=10;yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);xn00,fval,exitflag,output=fminbnd(yx,x11,x2)xn00=8.023562824723015fval=-3.568014059128578%最小值最小值exitflag=1output=iterations:9 funcCount:10 algorithm:golden section search,parabolic interpolation message:1x112 char y=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);（3）绘图观察最小值）绘图观察最小值xx=-10:pi/200:10;yxx=subs(y,x,xx);plot(xx,yxx)xlabel(x),grid on 12/6/2022第17页例4.1-8:f(x,y)=100(y-x2)2+(1-x)2在区间-5,5的极小值ff=(x)100*(x(2)-x(1).2).2+(1-x(1).2;%匿名函数 x0=-5,-2,2,5;-5,-2,2,5;%4个搜索起点 sx,sfval,sexit,soutput=fminsearch(ff,x0)其理论极小值点为其理论极小值点为x=1,y=1%sx给出一组使优化函数值非减的局部极小点给出一组使优化函数值非减的局部极小点 sx=0.99998 -0.68971 0.41507 8.0886 0.99997 -1.9168 4.9643 7.8004sfval=2.4112e-0102.4112e-010 5.7525e+002 2.2967e+003 3.3211e+005 format short e%取5位科学计数法 disp(ff(sx(:,1),ff(sx(:,2),ff(sx(:,3),ff(sx(:,4)用用x的的二元向量二元向量表示表示x,yx0=-5 -2 2 5 -5 -2 2 512/6/2022第18页4.1.5 常微分方程常微分方程Ordinary Differential Equation的数的数值值解解t,Y=ode45(odefun,tspan,y0)%4阶阶龙格库塔龙格库塔 数值数值法法odefun:待求解待求解一阶一阶微分方程组的函数文件句柄微分方程组的函数文件句柄tspan:自变量微分二元区间自变量微分二元区间t0,tfy0:一阶微分方程组的一阶微分方程组的(n*1)初值初值列向量列向量matlab为解常微分方程初值问题提供了一组配套齐全为解常微分方程初值问题提供了一组配套齐全,结构严结构严整的指令整的指令,包括包括:ode45,ode23,ode113,ode23t,ode15s,ode23s,ode23tb.在此只介绍最常用的在此只介绍最常用的ode45的基本使用方法的基本使用方法.ode45使用方法使用方法:t-二元区间二元区间的点系列的点系列Y-原函数在微分区间点系列上的函数值原函数在微分区间点系列上的函数值12/6/2022第19页例例4.1-9求解：求解：%解算微分方程解算微分方程tspan=0,30;y0=1;0;tt,yy=ode45(DyDt,tspan,y0);figure(1)plot(tt,yy(:,1)xlabel(t),title(x(t)解解:令令，上式写成，上式写成一阶一阶微分方程组形式微分方程组形式%画相平面图画相平面图(函数和其导数函数和其导数勾画的曲线称为勾画的曲线称为相轨迹相轨迹)figure(2)plot(yy(:,1),yy(:,2)xlabel(位移位移),ylabel(速速度度)function ydot=DyDt(t,y)mu=2;ydot=y(2);mu*(1-y(1)2)*y(2)-y(1);据以上方程组，编写据以上方程组，编写M函数文件函数文件DyDt.m12/6/2022第20页4.2矩阵和代数方程矩阵和代数方程4.2.1 矩阵运算和特征参数矩阵运算和特征参数4.2.2 矩阵的变换和特征值分解矩阵的变换和特征值分解4.2.3 线性方程的解线性方程的解4.2.4 一般代数方程的解一般代数方程的解12/6/2022第21页4.2.1 矩阵运算和特征参数矩阵运算和特征参数矩阵与标量之间的四则运算与数组运算相同矩阵与标量之间的四则运算与数组运算相同矩阵和矩阵之间的四则运算矩阵和矩阵之间的四则运算矩阵和矩阵之间的加减运算与数组运算相同矩阵和矩阵之间的加减运算与数组运算相同设设 A 是一个是一个 mn 矩阵，矩阵，B 是一个是一个 pq 矩阵，当矩阵，当 np 时，两个矩阵可以相乘，乘积为时，两个矩阵可以相乘，乘积为 mq 矩阵。矩阵矩阵。矩阵乘法不可逆。在乘法不可逆。在 MATLAB 中，矩阵乘法由中，矩阵乘法由“*”实实现。现。矩阵除法在实际中主要用于求解线性方程组矩阵除法在实际中主要用于求解线性方程组矩阵转置：矩阵转置：符号符号“”实现矩阵的转置操作。对于实数矩阵，实现矩阵的转置操作。对于实数矩阵，“”表示矩阵转置，对于复数矩阵，表示矩阵转置，对于复数矩阵，“”实现共轭转实现共轭转置。对于复数矩阵，如果想要实现非共轭转置，可以置。对于复数矩阵，如果想要实现非共轭转置，可以使用符号使用符号“.”。12/6/2022第22页format rat%有理格式显示有理格式显示A=magic(2)+j*pascal(2)A=1+1i 3+1i 4+1i 2+2iA1=A A1=1-1i 4-1i%共轭转置共轭转置 3-1i 2-2iA2=A.A2=1+1i 4+1i%非共轭转置，数组运算操作非共轭转置，数组运算操作 3+1i 2+2i例例4.2-2 矩阵和数组转置操作的差别矩阵和数组转置操作的差别12/6/2022第23页计算矩阵标量特征参数计算矩阵标量特征参数-秩秩,迹迹,行列式的指令行列式的指令rank(A)%求秩（求秩（Rank）det(A)%求行列式（求行列式（Determinant）trace(A)%求迹（求迹（Trace），即矩阵主对角元素的），即矩阵主对角元素的和和12/6/2022第24页 A=reshape(1:9,3,3);r=rank(A)%求秩求秩 d3=det(A)%非满秩矩阵的行列式一定为零非满秩矩阵的行列式一定为零 d2=det(A(1:2,1:2)%求子式的行列式求子式的行列式 t=trace(A)【例【例4.2-3】矩阵标量特征参数计算示例。】矩阵标量特征参数计算示例。A=1 4 7 2 5 8 3 6 9r=2d3=0d2=-3t=15 12/6/2022第25页【例【例4.2-4】矩阵标量特征参数的性质。】矩阵标量特征参数的性质。format short;rand(twister,0)A=rand(3,3);B=rand(3,3);C=rand(3,4);D=rand(4,3);tAB=trace(A*B)%任何符合矩阵乘法规则的两个矩阵乘积任何符合矩阵乘法规则的两个矩阵乘积tBA=trace(B*A)%的的“迹迹”不变。不变。同阶乘积同阶乘积“迹迹”不变不变tCD=trace(C*D)%两个两个“內维內维”相等矩阵的乘积相等矩阵的乘积“迹迹”不不变变tDC=trace(D*C)tAB=2.6030tBA=2.6030tCD=4.1191tDC=4.1191 dCD=det(C*D)dDC=det(D*C)dCD=0.0424 dDC=-2.6800e-018d_A_B=det(A)*det(B)dAB=det(A*B)dBA=det(B*A)同阶同阶矩阵乘积行列式等于各矩阵行列式之乘积矩阵乘积行列式等于各矩阵行列式之乘积d_A_B=0.0094dAB=0.0094dBA=0.0094 非同阶非同阶矩阵乘积行列矩阵乘积行列式式不等于不等于各矩阵行列各矩阵行列式之乘积式之乘积12/6/2022第26页udx=diff(X)%计计算向量算向量X的前向差分的前向差分uFX=gradient(F)%求一元求一元(函数函数)梯度梯度uFX,FY=gradient(F)%求二元（函数）梯度求二元（函数）梯度4.1.1 近似数值极限及导数近似数值极限及导数S1=sum(X,1)%sum按列向求和得按列向求和得(1n)数组数组,=sum(X)(X多行多行)S2=sum(X,2)%sum按行向求和得按行向求和得(n1)数组数组Scs=cumsum(X)%沿沿X列向求列向求累计和累计和,仍是数组仍是数组,最后一行等于最后一行等于Sx4.1.2 数值求和与近似数值积分数值求和与近似数值积分St=trapz(t,X)或或 St=dt*trapz(X)%梯形法求積分梯形法求積分Sct=cumtrapz(t,X)或或 Sct=dt*cumtrapz(X)%梯形法沿列向求梯形法沿列向求X关于关于x的的累计积分累计积分,最后一个值等于最后一个值等于StX=1 2;3 4 S1=4 6 S2=3 7sum(S1)=10sum(S1,1)=4 64 6 Review12/6/2022第27页4.1.3 计算精度可控的数值积分计算精度可控的数值积分一重积分一重积分(quadrature精度可控精度可控)：S1=quad(fun,a,b,tol)%自适应自适应Simpson法法S2=quadl(fun,a,b,tol)%自适应罗巴托自适应罗巴托 Lobatlo法法二重积分二重积分(精度可控精度可控)：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重积分三重积分(精度可控精度可控)：S4=triplequad(fun,xmin,xmax,ymin,ymax,zmix,zmax,tol)fun:可为字符串可为字符串,内联对象内联对象,匿名函数匿名函数,函数句柄函数句柄,注意乘注意乘除法运算一定加点除法运算一定加点.用数组运算用数组运算 tol:默认积分的绝对精度为默认积分的绝对精度为10-6Review12/6/2022第28页4.1.4 函数极函数极值值的数的数值值求解求解 x,fval,exitflag,output=fminbnd(fun,x1,x2,options)%一元函数一元函数在区间在区间bound(x1,x2)中极小值中极小值 x,fval:极值点坐标和对应目标函数极值极值点坐标和对应目标函数极值x,fval,exitflag,output=fminsearch(fun,x0,options)%单纯形法求搜索起点单纯形法求搜索起点x0附近附近多元函数多元函数极值点极值点%x每列代表一个候选极值点，各列按每列代表一个候选极值点，各列按目标函数极小值递增顺序目标函数极小值递增顺序%x(:,1)对应的目标函数极小值点由对应的目标函数极小值点由fval给出给出 fun:字符串字符串,内联函数内联函数,匿名函数匿名函数,函数句柄函数句柄,注意乘除法注意乘除法运算一定加点运算一定加点.用数组运算用数组运算 options:配置优化参数配置优化参数,可略可略 exitflag:给出大于给出大于0的数的数,则成功搜索到极值点则成功搜索到极值点 output:给出具体的优化算法和迭代次数给出具体的优化算法和迭代次数Review12/6/2022第29页rank(A)%求秩（求秩（Rank）det(A)%求行列式（求行列式（Determinant）trace(A)%求迹求迹(Trace)，即矩阵主对角元素的和，即矩阵主对角元素的和Reviewt,Y=ode45(odefun,tspan,y0)%4阶阶龙格库塔龙格库塔 数值数值法法odefun:待求解待求解一阶一阶微分方程组的函数文件句柄微分方程组的函数文件句柄tspan:自变量微分二元区间自变量微分二元区间t0,tfy0:一阶微分方程组的一阶微分方程组的(n*1)初值初值列向量列向量ode45使用方法使用方法:t-二元区间二元区间的点系列的点系列Y-原函数在微分区间点系列上的函数值原函数在微分区间点系列上的函数值4.1.5 常微分方程常微分方程Ordinary Differential Equation的数的数值值解解4.2.1 矩阵运算和特征参数矩阵运算和特征参数12/6/2022第30页4.2.2 矩阵的变换和特征值分解矩阵的变换和特征值分解R,ci=rref(A)%借助初等借助初等变换变换将将A缩减行变成行阶梯矩阵缩减行变成行阶梯矩阵R。B=orth(A)可得到矩阵可得到矩阵A的正交基，的正交基，B的列与的列与A的列可张成相同的列可张成相同的空间，而且的空间，而且B的列是正交的，因此的列是正交的，因此B*B=eye(rank(A)，B的的列数正好是列数正好是A的秩。的秩。X=null(A)%A矩阵零空间的全部正交基，满足矩阵零空间的全部正交基，满足AX=0,X*X=I。B=orth(A)V,D=eig(A)%A矩阵的特征值、特征向量分解，使矩阵的特征值、特征向量分解，使AV=VD。ci 是行数是行数组组,其元素表示其元素表示A中中线线性独立列的序号。性独立列的序号。length(ci)=rank(A)12/6/2022第31页【例【例4.2-5】行阶梯阵简化指令】行阶梯阵简化指令rref计算结果的含义。计算结果的含义。A=magic(4)R,ci=rref(A)%行阶梯分行阶梯分解解A=16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1R=1 0 0 1 0 1 0 3 0 0 1 -3 0 0 0 0ci=1 2 3 r_A=length(ci)%=rank(A)r_A=3 aa=A(:,1:3)*R(1:3,4)%A前三列线形组合前三列线形组合err=norm(A(:,4)-aa)aa=13 8 12 1err=0 12/6/2022第32页A=reshape(1:15,5,3)X=null(A)S=A*XS=1.0e-014*0 -0.1776 -0.2665 -0.3553 -0.5329A=1 6 11 2 7 12 3 8 13 4 9 14 5 10 15n=3l=1Rank_A=2X=0.4082 -0.8165 0.4082ans=1 例例4.2-6 矩矩阵零空零空间及含及含义。设设 是矩阵是矩阵 的零空间的零空间,即即n=size(A,2)l=size(X,2)Rank_A=rank(A)n-l=rank(A)12/6/2022第33页【例【例4.2-7】简单实阵简单实阵的特征分解。的特征分解。eig，cdf2rdf Complex Diagonal Form to Real-block Diagonal Form%（1）A=1,-3;2,2/3V,D=eig(A)A=1.0000 -3.0000 2.0000 0.6667V=0.7746 0.7746 0.0430-0.6310i 0.0430+0.6310iD=0.8333+2.4438i 0 0 0.8333-2.4438i%（2）VR,DR=cdf2rdf(V,D)VR=0.7746 0 0.0430 -0.6310DR=0.8333 2.4438 -2.4438 0.8333%（3）A1=V*D/V%因因计算算误差有很小虚部差有很小虚部A1_1=real(A1)%去除虚部后去除虚部后等于等于AA2=VR*DR/VRerr1=norm(A-A1,fro)err2=norm(A-A2,fro)A1=1.0000+0.0000i -3.0000 2.0000-0.0000i 0.6667 A1_1=1.0000 -3.0000 2.0000 0.6667A2=1.0000 -3.0000 2.0000 0.6667err1=6.7532e-016err2=4.4409e-016 Frobenius矩阵矩阵范数范数 实数块实数块对角阵对角阵实阵实阵矩阵矩阵Frobenius范数范数,类似向量类似向量2范数，不同于范数，不同于矩阵矩阵2范数范数12/6/2022第34页4.2.2 线性方程的解线性方程的解对于方程组对于方程组Amnx=b(m个方程，个方程，n个未知数个未知数)，当向量当向量b在矩阵在矩阵A列向量所张空间中列向量所张空间中,rank(A,b)=rank(A)=r a)若若n=r，则解唯一。，则解唯一。b)若若n r，则解不唯一。，则解不唯一。当向量当向量b不在矩阵不在矩阵A列向量所张空间中列向量所张空间中,则无准确解但则无准确解但存在最小二乘解。存在最小二乘解。1.线性方程解的一般讨论线性方程解的一般讨论matlab定义的左除运算可以很方便地解上述方程组：定义的左除运算可以很方便地解上述方程组：x=Ab (x=A-1b只适用于只适用于A非奇异时非奇异时)2.除法运算解除法运算解方程解方程解 当当m=n时时,“恰定恰定”方方程程 当当mn时时,“超定超定”方方程程 当当mr,解不唯一解不唯一%（3）求特解和通解，并对由之构成的全解进行验算）求特解和通解，并对由之构成的全解进行验算xs=Ab;xg=null(A);%xg是齐次方程是齐次方程Ax=0的解的解c=rand(1);ba=A*(xs+c*xg)%ba是是A乘乘“一个随机的全解一个随机的全解”Warning:Rank deficient,rank=2,tol=1.8757e-014.ba=13.0000 14.0000 15.0000 16.0000 例例4.2-8:求方程求方程 的的解解norm(ba-b)ans=1.1374e-014 12/6/2022第36页例例4.2-9:“逆阵法逆阵法”和和”左除法左除法”解恰定方程的性能对解恰定方程的性能对比。比。测试阵测试阵产生指定异常值和特殊带宽的随机阵产生指定异常值和特殊带宽的随机阵3.矩阵的逆矩阵的逆A_1=inv(A)%求非奇异阵求非奇异阵A的逆的逆randn(state,0);A=gallery(randsvd,300,2e13);%产生条件数产生条件数 为为2e13的的300阶随机矩阵阶随机矩阵;x=ones(300,1);%定义真解定义真解b=A*x;cond(A)%验算矩阵条件数验算矩阵条件数,结果结果a1.9978e+013,值越大值越大,阵越病态阵越病态,%求逆法求逆法tic;xi=inv(A)*b;ti=toceri=norm(x-xi)rei=norm(A*xi-b)/norm(b)%左除法左除法tic;xd=Ab;td=tocerd=norm(x-xd)red=norm(A*xd-b)/norm(b)ti=0.0185eri=0.0883rei=0.0051td=0.0066erd=0.0298red=8.7810e-015矩阵计算对于误差越敏感矩阵计算对于误差越敏感,数值稳定性数值稳定性差差 12/6/2022第37页求解任意函数求解任意函数f(x)=0(可能无解可能无解,单解单解,多解多解)的步骤的步骤:(1)作图获取初步近似解作图获取初步近似解 观察观察f(x)与横轴的交点坐标，用与横轴的交点坐标，用zoom放大，用放大，用ginput得较精确些的交点坐标值。得较精确些的交点坐标值。4.3 一般代数方程的解一般代数方程的解 (2)用用”泛函泛函”指令求精确解指令求精确解 x,favl=fzero(fun,x0)%一元函数求零点一元函数求零点 x,fval=fsolve(fun,x0)%解非解非线性方程性方程组 fun:内内联对象象,匿名函数匿名函数,函数句柄函数句柄,字符串字符串;被解函数自被解函数自变量一般用量一般用x,注意乘除法运算一定加点注意乘除法运算一定加点.用数用数组运算运算 x0:零点初始猜零点初始猜测值.x0为标量量时取与之最靠近的零点取与之最靠近的零点;x0取取a,b时,在此区在此区间内内寻找一个零点找一个零点.x:所求零点的自所求零点的自变量量值 fval:函数函数值12/6/2022第38页例4.2-10:求f(t)=(sin2t)e-0.1t0.5|t|的零点y_C=inline(sin(t).2.*exp(-0.1*t)-0.5*abs(t),t);t=-10:0.01:10;Y=y_C(t);clf,plot(t,Y,r);hold onplot(t,zeros(size(t),k);xlabel(t);ylabel(y(t)t4=0.5993y4=0tt=-2.0039 -0.5184 -0.0042 0.6052 1.6717 t3=-0.5198y3=5.5511e-017zoom ontt,yy=ginput(5),zoom offt4,y4=fzero(y_C,tt(4)t3,y3=fzero(y_C,tt(3)%t=0处没穿越横没穿越横轴12/6/2022第39页4.3 概率分布和统计分析概率分布和统计分析4.3.1 概率函数、分布函数、逆分布函概率函数、分布函数、逆分布函数和随机数的发生数和随机数的发生4.3.2 随机数发生器和随机数发生器和 统计分析指令统计分析指令12/6/2022第40页4.3 概率分布和统计分析概率分布和统计分析1.二项分布（二项分布（Binomial distribution）4.3.1 概率函数、分布函数、逆分布函数和随机数的发生概率函数、分布函数、逆分布函数和随机数的发生pk=binopdf(k,N,p)次的概率次的概率k发发生生A事件事件Fk=binocdf(k,N,p)的概率的概率k发发生次数不大于生次数不大于A事件事件R=binornd(N,p)产生符合二项分布产生符合二项分布B(N,p)的的(mn)随机随机数组（数组（元素值为元素值为事件事件A可能发生的次数可能发生的次数）。）。Binomial Cumulative Distribution FunctionBinomial Probability Density Function12/6/2022第41页N=100;p=0.5;k=0:N;pdf=binopdf(k,N,p);cdf=binocdf(k,N,p);h=plotyy(k,pdf,k,cdf);