《第4MATLAB的数值计算ppt课件》由会员分享，可在线阅读，更多相关《第4MATLAB的数值计算ppt课件（62页珍藏版）》请在金锄头文库上搜索。
1、第第4章章 MATLAB 的数值计算的数值计算 4.1 数值微积分数值微积分4.2 矩阵和代数方程矩阵和代数方程4.3 概率分布和统计分析概率分布和统计分析4.4 多项式运算和卷积多项式运算和卷积4.1数值微积分数值微积分4.1.1 近似数值极限及导数近似数值极限及导数4.1.2 数值求和与近似数值积分数值求和与近似数值积分4.1.3 计算精度可控的数值积分计算精度可控的数值积分4.1.4 函数极值的数值求解函数极值的数值求解4.1.5 常微分方程的数值解常微分方程的数值解u在在MATLAB数值计算中，既没有专门的求极限指令，数值计算中，既没有专门的求极限指令，也没有专门的求导指令。但也没有专
2、门的求导指令。但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 近似数值极限及导数近似数值极限及导数udx=diff(
3、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 , u FX(end) = F(end)-F(end-1) ; FX长度与长度与F的长度一样的长度一样u当当F是矩阵时，是矩阵时， FX, FY是与是与F同样大小的矩阵。同样大小的矩阵。 FX的的每行对应每行对应F相应行元素间的梯度相应行元素间
4、的梯度 ; FY的每列对应的每列对应F相应相应列元素间的梯度列元素间的梯度 ; 4.1.1 近似数值极限及导数近似数值极限及导数数值极限和导数的运用应非常谨慎数值极限和导数的运用应非常谨慎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; %可得到正确结果可得到正确结果数值极限和导数的运用应非常谨慎数值极限和导数的运用应非常谨慎%
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)的数的数值非常接近，高位有效数字完
6、全一非常接近，高位有效数字完全一样， df =f(t+dt)-f(t) 呵斥高位有效数字消逝，精度呵斥高位有效数字消逝，精度急急剧变差。差。数值极限和导数的运用应非常谨慎数值极限和导数的运用应非常谨慎%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;d=pi/100; t=0:d:2*
7、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,Loc
8、ation,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结果大致一样结果大
9、致一样细节上细节上, diff和和gradient数值有差别，数值有差别， diff没有给出最后一点导数没有给出最后一点导数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*cumtra
10、pz(X) %梯形法沿列向求梯形法沿列向求X关于关于x的累计积分的累计积分,最后一最后一个值等于个值等于StS=dt*sum(X), S=sum(t,X) %近似矩形法求近似矩形法求积分积分Scs=cumsum(t,X) =dt*cumsum(X) 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
11、求得积分求得积分)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.30134.1.3 计算精度可控的数值积分计算精度可控的数值积分一重积分一重积分(quadrature精度可控精度可控)：S1=quad(fun,a,b,tol) %自顺应自顺
12、应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-61syms
13、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 3fx=exp(-x.2); %一定用数组乘一定用数组乘Ic=quad(fx,0,1,1e-8) Ic = 0.746824132854452 %实践精度控制到实践精度控
14、制到1e-101符号计算法符号计算法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=
15、dblquad(inline(x.y),0,1,1,2) %内联函数表示被积函数内联函数表示被积函数% 一定为数组乘一定为数组乘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附近多元函数极值点附近多元函数极值点
16、 % x每列代表一个候选极值点，各列按目的函数每列代表一个候选极值点，各列按目的函数极小值递增顺序极小值递增顺序 %x(:,1)对应的目的函数极小值点由对应的目的函数极小值点由fval给出给出 fun: 字符串字符串,内联函数内联函数,匿名函数匿名函数,函数句柄函数句柄,留留意乘除法运算一定加点意乘除法运算一定加点.用数组运算用数组运算 options: 配置优化参数配置优化参数,可略可略 exitflag: 给出大于给出大于0的数的数,那么胜利搜索到极值点那么胜利搜索到极值点 output: 给出详细的优化算法和迭代次数给出详细的优化算法和迭代次数【例【例4.1-7 】知】知 ，在，在-10
17、x10区区间，求函数的最小，求函数的最小值。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
18、 = 2.2958874039497802890014385492622*10(-41)y_xs0 = -0.001263317776974196724544154344118 无法判别能否最小值无法判别能否最小值【例【例4.1-7 】知】知 ，在，在-10x10区区间，求函数的最小，求函数的最小值。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.514797840754235fva
19、l = %比比“导数为零法求得的极值更小导数为零法求得的极值更小, 更能够是最小更能够是最小值值 -0.499312445280039exitflag = 1output = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation message: 1x112 char 【例【例4.1-7 】知】知 ，在，在-10x10区区间，求函数的最小，求函数的最小值。4据图形察看，重设据图形察看，重设fminbnd的搜索区间的搜索区间x11=6;x2=10;yx=(x)(sin(x)2*
20、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)
21、;3绘图察看最小值绘图察看最小值xx=-10:pi/200:10;yxx=subs(y,x,xx);plot(xx,yxx)xlabel(x),grid on 例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给出一组使优化函数值非减的部分极小点给出一组使优化函数值非减的部分极小
22、点 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 54.1.5 常微分方程常微分方程Ordinary Differential Equation的的
23、数数值解解t,Y=ode45(odefun,tspan,y0) % 4阶龙格库塔阶龙格库塔 数值法数值法odefun: 待求解一阶微分方程组的函数文件句柄待求解一阶微分方程组的函数文件句柄tspan: 自变量微分二元区间自变量微分二元区间t0, tfy0: 一阶微分方程组的一阶微分方程组的(n*1)初值列向量初值列向量matlab为解常微分方程初值问题提供了一组配套齐全为解常微分方程初值问题提供了一组配套齐全,构造严构造严整的指令整的指令, 包括包括: ode45, ode23, ode113, ode23t, ode15s, ode23s, ode23tb.在此只引见最常用的在此只引见最常用
24、的ode45的根本运用方法的根本运用方法.ode45运用方法运用方法:t-二元区间的点系列二元区间的点系列Y-原函数在微分区间点系列上的函数值原函数在微分区间点系列上的函数值例例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(:
25、,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.m4.2矩阵和代数方程矩阵和代数方程4.2.1 矩阵运算和特征参数矩阵运算和特征参数4.2.2 矩阵的变换和特征值分解矩阵的变换和特征值分解4.2.3 线性方程的解线性方程的解4.2.4 普通代数方程的解普通代数方程的解4.2.1 矩阵运算和特征参数矩阵运算和特征参数矩阵与标量之间的四那么运算与数组运算一样矩阵与标量之间的四那么运算与数
26、组运算一样矩阵和矩阵之间的四那么运算矩阵和矩阵之间的四那么运算矩阵和矩阵之间的加减运算与数组运算一样矩阵和矩阵之间的加减运算与数组运算一样设设 A 是一个是一个 mn 矩阵，矩阵，B 是一个是一个 pq 矩阵，当矩阵，当 np 时，两个矩阵可以相乘，乘积为时，两个矩阵可以相乘，乘积为 mq 矩阵。矩阵乘法矩阵。矩阵乘法不可逆。在不可逆。在 MATLAB 中，矩阵乘法由中，矩阵乘法由“*实现。实现。矩阵除法在实践中主要用于求解线性方程组矩阵除法在实践中主要用于求解线性方程组矩阵转置：矩阵转置：符号符号“ 实现矩阵的转置操作。对于实数矩阵，实现矩阵的转置操作。对于实数矩阵， “ 表示矩阵转置，对于
27、复数矩阵，表示矩阵转置，对于复数矩阵，“ 实现共轭转实现共轭转置。对于复数矩阵，假设想要实现非共轭转置，可以置。对于复数矩阵，假设想要实现非共轭转置，可以运用符号运用符号“ . 。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 矩阵和数组转置操作的差别矩阵和
28、数组转置操作的差别计算矩阵标量特征参数计算矩阵标量特征参数-秩秩,迹迹,行列式的指令行列式的指令rank(A) %求秩求秩Rankdet(A) %求行列式求行列式Determinanttrace(A) %求迹求迹Trace，即矩阵主对角元素，即矩阵主对角元素的和的和 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
29、 5 8 3 6 9r = 2d3 =0d2 = -3t = 15 【例【例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
30、*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 非同阶矩阵乘积行列非同阶矩阵乘积行列式不等于各矩阵行列式不等于各矩阵行列式之乘积式之乘积udx=diff(X) %计算向量算向量X的前向
31、差分的前向差分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*tra
32、pz(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 6Review4.1.3 计算精度可控的数值积分计算精度可控的数值积分一重积分一重积分(quadrature精度可控精度可控)：S1=quad(fun,a,b,tol) %自顺应自顺应Simpson法法S2=quadl(fun,a,b,tol) %自顺应罗巴托自顺应罗巴托 Lob
33、atlo法法二重积分二重积分(精度可控精度可控) ：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重积分三重积分(精度可控精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可为字符串可为字符串,内联对象内联对象,匿名函数匿名函数,函数句柄函数句柄,留意乘留意乘除法运算一定加点除法运算一定加点.用数组运算用数组运算 tol: 默许积分的绝对精度为默许积分的绝对精度为10-6Review4.1.4 函数极值的数值求解函数极值的数值求解 x,fval,exitflag,output=f
34、minbnd(fun,x1,x2,options) %一元函数在区间一元函数在区间bound(x1,x2)中极小值中极小值 x,fval: 极值点坐标和对应目的函数极值极值点坐标和对应目的函数极值x,fval,exitflag,output=fminsearch(fun,x0,options) %单纯形法求搜索起点单纯形法求搜索起点x0附近多元函数极值点附近多元函数极值点 % x每列代表一个候选极值点，各列按目的函数每列代表一个候选极值点，各列按目的函数极小值递增顺序极小值递增顺序 %x(:,1)对应的目的函数极小值点由对应的目的函数极小值点由fval给出给出 fun: 字符串字符串,内联函数
35、内联函数,匿名函数匿名函数,函数句柄函数句柄, 留留意乘除法运算一定加点意乘除法运算一定加点.用数组运算用数组运算 options: 配置优化参数配置优化参数,可略可略 exitflag: 给出大于给出大于0的数的数,那么胜利搜索到极值点那么胜利搜索到极值点 output: 给出详细的优化算法和迭代次数给出详细的优化算法和迭代次数Reviewrank(A) %求秩求秩Rankdet(A) %求行列式求行列式Determinanttrace(A) %求迹求迹(Trace)，即矩阵主对角元素的和，即矩阵主对角元素的和Reviewt,Y=ode45(odefun,tspan,y0) % 4阶龙格库塔
36、阶龙格库塔 数值法数值法odefun: 待求解一阶微分方程组的函数文件句柄待求解一阶微分方程组的函数文件句柄tspan: 自变量微分二元区间自变量微分二元区间t0, tfy0: 一阶微分方程组的一阶微分方程组的(n*1)初值列向量初值列向量ode45运用方法运用方法:t-二元区间的点系列二元区间的点系列Y-原函数在微分区间点系列上的函数值原函数在微分区间点系列上的函数值4.1.5 常微分方程常微分方程Ordinary Differential Equation的的数数值解解4.2.1 矩阵运算和特征参数矩阵运算和特征参数4.2.2 矩阵的变换和特征值分解矩阵的变换和特征值分解R, ci=rre
37、f(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中中线性独立列的序号。
38、性独立列的序号。length(ci)=rank(A)【例【例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 1
39、err = 0 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)【例【例4.2-7】简单实阵的特征分解。的特征分解。eig，cdf2rdf
40、Complex Diagonal Form to Real-block Diagonal Form%1A=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 %2VR,DR=cdf2rdf(V,D) VR =0.7746 0 0.0430 -0.6310DR =0.8333 2.4438 -2.4438 0.8333 %3A1=V*D/V %因因计算算误差有很差
41、有很小虚部小虚部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范数范数,类似向量类似
42、向量2范数，不同于范数，不同于矩阵矩阵2范数范数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定义的左除运算可以很方便地解上述方程组
43、：定义的左除运算可以很方便地解上述方程组： 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.
44、8757e-014.ba =13.0000 14.0000 15.0000 16.0000 例例4.2-8: 求方程求方程 的解的解norm(ba-b) ans = 1.4e-014 例例4.2-9: “逆阵法和左除法解恰定方程的性能对逆阵法和左除法解恰定方程的性能对比。比。测试阵测试阵产生指定异常值和特殊带宽的随机阵产生指定异常值和特殊带宽的随机阵3. 矩阵的逆矩阵的逆A_1=inv(A) %求非奇特阵求非奇特阵A的逆的逆randn(state,0); A=gallery(randsvd,300,2e13); %产生条件数产生条件数 为为2e13的的300阶随机矩阵阶随机矩阵;x=ones(
45、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矩阵计算对于误差越敏感矩阵计算
46、对于误差越敏感, 数值稳定性数值稳定性差差 求解恣意函数求解恣意函数f(x)=0(能够无解能够无解,单解单解,多解多解)的步骤的步骤:作图获取初步近似解作图获取初步近似解 察看察看f(x)与横轴的交点坐标，用与横轴的交点坐标，用zoom放大，用放大，用ginput得较准确些的交点坐标值。得较准确些的交点坐标值。4.3 普通代数方程的解普通代数方程的解 (2)用泛函指令求准确解用泛函指令求准确解(3) x,favl = fzero(fun,x0) %一元函数求零点一元函数求零点(4) x,fval = fsolve(fun,x0) %解非解非线性方程性方程组(5) fun: 内内联对象象,匿名函
47、数匿名函数,函数句柄函数句柄,字符串字符串;被解函数被解函数自自变量普通用量普通用x, 留意乘除法运算一定加点留意乘除法运算一定加点.用数用数组运运算算(6) x0: 零点初始猜零点初始猜测值. x0为标量量时取与之最接近的零取与之最接近的零点点;(7) x0取取a,b时,在此区在此区间内内寻觅一个零点一个零点.(8) x: 所求零点的自所求零点的自变量量值(9) fval: 函数函数值例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(
48、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处没穿越没穿越横横轴4.3 概率分布和统计分析概率分布和统计分析4.3.1 概率函数、分布函数、逆分布概率函数、分布函
49、数、逆分布函数和随机数的发生函数和随机数的发生4.3.2 随机数发生器和随机数发生器和 统计分析指统计分析指令令4.3 概率分布和统计分析概率分布和统计分析1. 二项分布二项分布Binomial distribution 4.3.1 概率函数、分布函数、逆分布函数和随机数的发生概率函数、分布函数、逆分布函数和随机数的发生次的概率次的概率k发生生A事件事件pk=binopdf(k,N,p) 的的k发生次数不大于生次数不大于A事件事件Fk=binocdf(k,N,p) 概率概率R=binornd(N,p) 产生符合二项分布产生符合二项分布B(N,p)的的(mn) 随机随机数组元素值为事件数组元素值
50、为事件A能够发生的次数能够发生的次数。Binomial Cumulative Distribution FunctionBinomial Probability Density FunctionN=100;p=0.5;k=0:N;pdf=binopdf(k,N,p); cdf=binocdf(k,N,p);h=plotyy(k,pdf,k,cdf);set(get(h(1),Children),Color,b,Marker,.,MarkerSize,13)set(get(h(1),Ylabel),String,pdf)set(h(2),Ycolor,1,0,0)set(get(h(2),Chi
51、ldren),Color,r,Marker,+,MarkerSize,4)set(get(h(2),Ylabel),String,cdf)xlabel(k) grid on 例例4.3-1. 画出画出N=100，p=0.5情况下的二项分布概率特性情况下的二项分布概率特性曲线。曲线。4.3.1 概率函数、分布函数、逆分布函数和随机数的发生概率函数、分布函数、逆分布函数和随机数的发生px=normpdf(x,Mu,Sigma) Fx=normcdf(x,Mu,Sigma)R=normrnd(Mu,Sigma,m,n)2. 正态分布正态分布Normal distribution 的概率密度。的概率密
52、度。x分布的随机分布的随机变量取量取值 服从服从的概率分布。的概率分布。x分布的随机分布的随机变量取量取值不大于不大于 服从服从产生元素服从产生元素服从 分布的分布的(mn) 随机数组。随机数组。, 是正是正态分布的数学期望分布的数学期望Mu可取任何可取任何实数，数，x是正是正态分布的均方差。分布的均方差。Sigma 例例4.3-2. 正态分布规范正态分布规范差的几何表示。差的几何表示。mu=3;sigma=0.5;x=mu+sigma*-3:-1,1:3yf=normcdf(x,mu,sigma)x = 1.5 2.0 2.5 3.5 4.0 4.5yf = 0.0013 0.0228 0.
53、1587 0.8413 0.9772 0.9987P=yf(4)-yf(3),yf(5)-yf(2),yf(6)-yf(1) P = 0.6827 0.9545 0.9973xd=1:0.1:5; yd=normpdf(xd,mu,sigma);例例4.3-2. 正态分布规范正态分布规范差的几何表示。差的几何表示。clffor k=1:3xx=x(4-k):sigma/10:x(3+k);yy=normpdf(xx,mu,sigma);subplot(3,1,k),plot(xd,yd,b);hold onfill(x(4-k),xx,x(3+k),0,yy,0,g);hold offif k
54、 Open Unable to read MAT file C:MATLAB6p5p1work prob_data401.matFile may be corrupt 由于由于mat文件格式的变化不兼容，低版本无法翻开高文件格式的变化不兼容，低版本无法翻开高版本的文件很常见。详细缘由如下：从版本的文件很常见。详细缘由如下：从7.0起新的起新的MAT格式中，数字量先紧缩再保管，字符串也是以格式中，数字量先紧缩再保管，字符串也是以Unicode编码保管的，默许情况下保管为新格式，所编码保管的，默许情况下保管为新格式，所以只能在以只能在7.X中翻开。为兼容他可以重新运用中翻开。为兼容他可以重新运用
55、save filename (文件名文件名) -v6 保管为兼容方式，就行保管为兼容方式，就行了。了。 Matlab安装乱码安装乱码-处理方案处理方案处理方案：单击处理方案：单击“File“Preferences“Fonts“Desktop code font栏中选栏中选“微软雅黑微软雅黑, Plain, 10“Custom fonts每一项包括每一项包括“command window都点选都点选“Desktop code即可即可.本人亦可选择本人中意的字体与本人亦可选择本人中意的字体与custom fonts方案，方案， 系统：系统：Win7Ultimate x86安装安装Matlab版本：
56、版本：Matlab R2019b问题问题：在安装好：在安装好R2019b之后，启动试用，结果翻开之后，启动试用，结果翻开发现发现Command Window全是乱码，在里面输入命令也全是乱码，在里面输入命令也全是一堆乱符，如下所示全是一堆乱符，如下所示.习题习题4 (Page186) 重点重点2, 5, 6, 9,10,12,17,18 题题留意留意4, matlab6.5符号计算计算结果大符号计算计算结果大1, s = 2.0878494, 应为应为s = 1.0878494, ，不做不做 7，8，11, 13(涉及涉及contour在后面在后面), 14,15,1611题为超定方程，独立方程个数大于独立的未知参数的题为超定方程，独立方程个数大于独立的未知参数的个数的方程，在个数的方程，在matlab里面有三种方法求解，一是用伪里面有三种方法求解，一是用伪逆法求解，逆法求解，x=pinv(A)*b,二是用左除法求解，二是用左除法求解，x=Ab,三是三是用最小二乘法求解用最小二乘法求解,x=lsqnonneg(A,b)
