第4MATLAB的數(shù)值計(jì)算-ppt課件

1、第第4章章 MATLAB 的數(shù)值計(jì)算的數(shù)值計(jì)算 4.1 數(shù)值微積分?jǐn)?shù)值微積分4.2 矩陣和代數(shù)方程矩陣和代數(shù)方程4.3 概率分布和統(tǒng)計(jì)分析概率分布和統(tǒng)計(jì)分析4.4 多項(xiàng)式運(yùn)算和卷積多項(xiàng)式運(yùn)算和卷積4.1數(shù)值微積分?jǐn)?shù)值微積分4.1.1 近似數(shù)值極限及導(dǎo)數(shù)近似數(shù)值極限及導(dǎo)數(shù)4.1.2 數(shù)值求和與近似數(shù)值積分?jǐn)?shù)值求和與近似數(shù)值積分4.1.3 計(jì)算精度可控的數(shù)值積分計(jì)算精度可控的數(shù)值積分4.1.4 函數(shù)極值的數(shù)值求解函數(shù)極值的數(shù)值求解4.1.5 常微分方程的數(shù)值解常微分方程的數(shù)值解u在在MATLAB數(shù)值計(jì)算中，既沒有專門的求極限指令，數(shù)值計(jì)算中，既沒有專門的求極限指令，也沒有專門的求導(dǎo)指令。但也沒有專

2、門的求導(dǎo)指令。但MATLAB提供了與提供了與“求導(dǎo)求導(dǎo)概念有關(guān)的概念有關(guān)的“求差分指令。求差分指令。udx=diff(X) %計(jì)算向量計(jì)算向量X的前向差分的前向差分uFX=gradient(F) %求一元求一元(函數(shù)函數(shù))梯度梯度uFX, FY =gradient(F) %求二元函數(shù)梯度求二元函數(shù)梯度u對對diff而言，當(dāng)而言，當(dāng)X是向量時(shí)，是向量時(shí)，dx= X(2:n)-X(1:n-1) ;當(dāng)當(dāng)X是矩陣時(shí)，是矩陣時(shí)，dx= X(2:n, :)-X(1:n-1, :) 。 dx的長度比的長度比x的的長度少長度少1個(gè)元素。個(gè)元素。4.1.1 近似數(shù)值極限及導(dǎo)數(shù)近似數(shù)值極限及導(dǎo)數(shù)udx=diff(

3、X) %計(jì)算向量計(jì)算向量X的前向差分的前向差分uFX=gradient(F) %求一元求一元(函數(shù)函數(shù))梯度梯度uFX, FY =gradient(F) %求二元函數(shù)梯度求二元函數(shù)梯度u對對gradient而言，當(dāng)而言，當(dāng)F是向量時(shí)，是向量時(shí)，F(xiàn)X(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當(dāng)當(dāng)F是矩陣時(shí)，是矩陣時(shí)， FX, FY是與是與F同樣大小的矩陣。同樣大小的矩陣。 FX的的每行對應(yīng)每行對應(yīng)F相應(yīng)行元素間的梯度相應(yīng)行元素間

4、的梯度 ; FY的每列對應(yīng)的每列對應(yīng)F相應(yīng)相應(yīng)列元素間的梯度列元素間的梯度 ; 4.1.1 近似數(shù)值極限及導(dǎo)數(shù)近似數(shù)值極限及導(dǎo)數(shù)數(shù)值極限和導(dǎo)數(shù)的運(yùn)用應(yīng)非常謹(jǐn)慎數(shù)值極限和導(dǎo)數(shù)的運(yùn)用應(yīng)非常謹(jǐn)慎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; %可得到正確結(jié)果可得到正確結(jié)果數(shù)值極限和導(dǎo)數(shù)的運(yùn)用應(yīng)非常謹(jǐn)慎數(shù)值極限和導(dǎo)數(shù)的運(yùn)用應(yīng)非常謹(jǐn)慎%

5、1自變量的增量獲得過小自變量的增量獲得過小eps數(shù)量級數(shù)量級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) )sin(tx 2 , 0【例【例4.1-2】知】知, 求該函數(shù)在區(qū)間求該函數(shù)在區(qū)間 中的近似導(dǎo)函數(shù)。中的近似導(dǎo)函數(shù)。01234567-1.5-1-0.500.51t x(t)dx/dt數(shù)值導(dǎo)數(shù)受計(jì)算中有限精度困擾，當(dāng)增量數(shù)

6、值導(dǎo)數(shù)受計(jì)算中有限精度困擾，當(dāng)增量dt過小時(shí)，過小時(shí)，f(t+dt)與與f(t)的數(shù)值非常接近，高位有效數(shù)字完全一的數(shù)值非常接近，高位有效數(shù)字完全一樣，樣， df =f(t+dt)-f(t) 呵斥高位有效數(shù)字消逝，精度呵斥高位有效數(shù)字消逝，精度急劇變差。急劇變差。數(shù)值極限和導(dǎo)數(shù)的運(yùn)用應(yīng)非常謹(jǐn)慎數(shù)值極限和導(dǎo)數(shù)的運(yùn)用應(yīng)非常謹(jǐn)慎%2自變量的增量獲得適當(dāng)自變量的增量獲得適當(dāng)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) )sin(tx

7、2 , 0【例【例4.1-2】知】知, 求該函數(shù)在區(qū)間求該函數(shù)在區(qū)間 中的近似導(dǎo)函數(shù)。中的近似導(dǎo)函數(shù)。01234567-1-0.8-0.6-0.4-0.200.20.40.60.81t x(t)dx/dtd=pi/100;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計(jì)算該函計(jì)算該函數(shù)在區(qū)間數(shù)在區(qū)間 中的近似導(dǎo)函數(shù)。中的近似導(dǎo)函數(shù)。)sin( tx 2 , 0subplot(1,2,1); plot(t,x,b);hold onp

8、lot(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

9、-1),.k,MarkerSize,8)title(end-10, end)legend(dxdt_grad,dxdt_diff,Location,SouthEast)xlabel(t),box off; hold off 宏觀上宏觀上, diff和和gradient結(jié)果大致一樣結(jié)果大致一樣細(xì)節(jié)上細(xì)節(jié)上, diff和和gradient數(shù)值有差別，數(shù)值有差別， diff沒有給出最后一點(diǎn)導(dǎo)數(shù)沒有給出最后一點(diǎn)導(dǎo)數(shù)Sx=sum(X) % sum按列向求和得按列向求和得(1n)數(shù)組數(shù)組Scs=cumsum(X)%沿沿X列向求累計(jì)和列向求累計(jì)和, 仍是數(shù)組仍是數(shù)組, 第第(i, k)個(gè)元素是個(gè)元素是X數(shù)組

10、第數(shù)組第k列前列前i個(gè)元素的和。最后一行等于個(gè)元素的和。最后一行等于Sx4.1.2 數(shù)值求和與近似數(shù)值積分?jǐn)?shù)值求和與近似數(shù)值積分1( )( ,)mm niSx kXi k111221223132xxXxxxxSt=trapz(t,X) 或或 St=dt*trapz(X) %梯形法求積梯形法求積分分Sct=cumtrapz(t,X) 或或 Sct=dt*cumtrapz(X) %梯形法沿列向求梯形法沿列向求X關(guān)于關(guān)于x的累計(jì)積分的累計(jì)積分,最后一最后一個(gè)值等于個(gè)值等于St111211211222112131122232xxScsxxxxxxxxxx112112221121213112222232

11、00()/2()/2()/2 ()/2SctxxxxxxxxxxxxS=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); %梯形法積分梯形法積分2/0)()(dttyxs)sin(2 . 0ty【例【例4.1-4】求積分】求積分 , 其中其中disp(sum求得積分求得積分,blanks(3),

12、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.301300.20.40.60.811.21.41.61.800.511.54.1.3 計(jì)算精度可控的數(shù)值積分計(jì)算精度可控的數(shù)值積分一重積分一重積分(quadratu

13、re精度可控精度可控)：S1=quad(fun,a,b,tol) %自順應(yīng)自順應(yīng)Simpson法法S2=quadl(fun,a,b,tol) %自順應(yīng)羅巴托自順應(yīng)羅巴托 Lobatlo法法二重積分二重積分(精度可控精度可控) ：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重積分三重積分(精度可控精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可為字符串可為字符串,內(nèi)聯(lián)對象內(nèi)聯(lián)對象,匿名函數(shù)匿名函數(shù),函數(shù)句柄函數(shù)句柄,留意乘除留意乘除法運(yùn)算一定加點(diǎn)法運(yùn)算一定加點(diǎn).用數(shù)組運(yùn)算用數(shù)

14、組運(yùn)算 tol: 默許積分的絕對精度為默許積分的絕對精度為10-61syms xIsym=vpa(int(exp(-x2),x,0,1) % x.2數(shù)組乘方亦可數(shù)組乘方亦可Isym = 0.74682413281242702539946743613185 例例4.1-5：求：求 . 102dxeIx% 2 梯形法積分梯形法積分format longd=0.001;x=0:d:1;Itrapz=d*trapz(exp(-x.*x) % x.*x必需為數(shù)組乘必需為數(shù)組乘Itrapz = 0.746824071499185 3fx=exp(-x.2); %一定用數(shù)組乘一定用數(shù)組乘Ic=quad(fx

15、,0,1,1e-8) Ic = 0.746824132854452 %實(shí)踐精度控制到實(shí)踐精度控制到1e-101符號計(jì)算法符號計(jì)算法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：求：求 . 2110ysx dxdy 2數(shù)值積分法數(shù)值積分法format longs_n=dblquad(x,y)x.y,0,1,1,2) %匿名函數(shù)表示被積函數(shù)匿名函數(shù)表示被積函數(shù)s_n = 0.4054

16、66267243508 s_n=dblquad(x.y,0,1,1,2) %字符串表示被積函數(shù)字符串表示被積函數(shù)s_n=dblquad(inline(x.y),0,1,1,2) %內(nèi)聯(lián)函數(shù)表示被積函數(shù)內(nèi)聯(lián)函數(shù)表示被積函數(shù)% 一定為數(shù)組乘一定為數(shù)組乘4.1.4 函數(shù)極值的數(shù)值求解函數(shù)極值的數(shù)值求解 x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函數(shù)在區(qū)間一元函數(shù)在區(qū)間bound(x1,x2)中極小值中極小值 x,fval: 極值點(diǎn)坐標(biāo)和對應(yīng)目的函數(shù)極值極值點(diǎn)坐標(biāo)和對應(yīng)目的函數(shù)極值x,fval,exitflag,output=fminse

17、arch(fun,x0,options) %單純形法求搜索起點(diǎn)單純形法求搜索起點(diǎn)x0附近多元函數(shù)極值點(diǎn)附近多元函數(shù)極值點(diǎn) % x每列代表一個(gè)候選極值點(diǎn)，各列按目的函數(shù)每列代表一個(gè)候選極值點(diǎn)，各列按目的函數(shù)極小值遞增順序極小值遞增順序 %x(:,1)對應(yīng)的目的函數(shù)極小值點(diǎn)由對應(yīng)的目的函數(shù)極小值點(diǎn)由fval給出給出 fun: 字符串字符串,內(nèi)聯(lián)函數(shù)內(nèi)聯(lián)函數(shù),匿名函數(shù)匿名函數(shù),函數(shù)句柄函數(shù)句柄,留意留意乘除法運(yùn)算一定加點(diǎn)乘除法運(yùn)算一定加點(diǎn).用數(shù)組運(yùn)算用數(shù)組運(yùn)算 options: 配置優(yōu)化參數(shù)配置優(yōu)化參數(shù),可略可略 exitflag: 給出大于給出大于0的數(shù)的數(shù),那么勝利搜索到極值點(diǎn)那么勝利搜索到極值

18、點(diǎn) output: 給出詳細(xì)的優(yōu)化算法和迭代次數(shù)給出詳細(xì)的優(yōu)化算法和迭代次數(shù)【例【例4.1-7 】知】知 ，在，在-10 x10區(qū)間，求函數(shù)的最小值。區(qū)間，求函數(shù)的最小值。0.12sin0.5(0.1)sinxyexxx1用用“導(dǎo)數(shù)為零法求極值點(diǎn)導(dǎo)數(shù)為零法求極值點(diǎn)syms xy=sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);yd=diff(y,x); %求導(dǎo)函數(shù)求導(dǎo)函數(shù)xs0=solve(yd,x) %求導(dǎo)函數(shù)為零的根，即極值點(diǎn)求導(dǎo)函數(shù)為零的根，即極值點(diǎn)yd_xs0=vpa(subs(yd,x,xs0) %驗(yàn)證極值點(diǎn)處導(dǎo)驗(yàn)證極值點(diǎn)處導(dǎo)函數(shù)為零函數(shù)為零y_xs0=

19、vpa(subs(y,x,xs0) %求極值點(diǎn)處極值求極值點(diǎn)處極值xs0 = matrix(0.050838341410271656880659496266968)yd_xs0 = 2.2958874039497802890014385492622*10(-41)y_xs0 = -0.001263317776974196724544154344118 無法判別能否最小值無法判別能否最小值【例【例4.1-7 】知】知 ，在，在-10 x10區(qū)間，求函數(shù)的最小值。區(qū)間，求函數(shù)的最小值。0.12sin0.5(0.1)sinxyexxx2采用優(yōu)化算法求極小值采用優(yōu)化算法求極小值x1=-10;x2=10

20、;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 = %比比“導(dǎo)數(shù)為零法求得的極值更小導(dǎo)數(shù)為零法求得的極值更小, 更能夠是最小更能夠是最小值值 -0.499312445280039exitflag = 1output = iterations: 13 funcCount: 14 algorithm: golden section search, parabolic interpolation message:

21、1x112 char 【例【例4.1-7 】知】知 ，在，在-10 x10區(qū)間，求函數(shù)的最小值。區(qū)間，求函數(shù)的最小值。0.12sin0.5(0.1)sinxyexxx4據(jù)圖形察看，重設(shè)據(jù)圖形察看，重設(shè)fminbnd的搜索區(qū)間的搜索區(qū)間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 = it

22、erations: 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 -10-8-6-4-20246810-4-3-2-101234x例4.1-8: f(x,y)=100(y-x2)2+(1-x)2在區(qū)間-5,5的極小值f

23、f=(x)100*(x(2)-x(1).2).2+(1-x(1).2; %匿名函數(shù) x0=-5,-2,2,5;-5,-2,2,5; %4個(gè)搜索起點(diǎn) sx,sfval,sexit,soutput=fminsearch(ff,x0)其實(shí)際極小值點(diǎn)為其實(shí)際極小值點(diǎn)為x=1, y=1%sx給出一組使優(yōu)化函數(shù)值非減的部分極小點(diǎn)給出一組使優(yōu)化函數(shù)值非減的部分極小點(diǎn) 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

24、.3211e+005 format short e %取5位科學(xué)計(jì)數(shù)法 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的的數(shù)值解數(shù)值解t,Y=ode45(odefun,tspan,y0) % 4階龍格庫塔階龍格庫塔 數(shù)數(shù)值法值法odefun: 待求解一階微分方程組的函數(shù)文件句柄待求解一階微分方程組的函數(shù)文件句柄tspan: 自變量微分二元區(qū)間自變量微分二元區(qū)間t0,

25、tfy0: 一階微分方程組的一階微分方程組的(n*1)初值列向量初值列向量matlab為解常微分方程初值問題提供了一組配套齊全為解常微分方程初值問題提供了一組配套齊全,構(gòu)造嚴(yán)構(gòu)造嚴(yán)整的指令整的指令, 包括包括: ode45, ode23, ode113, ode23t, ode15s, ode23s, ode23tb.在此只引見最常用的在此只引見最常用的ode45的根本運(yùn)用方法的根本運(yùn)用方法.ode45運(yùn)用方法運(yùn)用方法:t-二元區(qū)間的點(diǎn)系列二元區(qū)間的點(diǎn)系列Y-原函數(shù)在微分區(qū)間點(diǎn)系列上的函數(shù)值原函數(shù)在微分區(qū)間點(diǎn)系列上的函數(shù)值例例4.1-9求解：求解：222(0)(1)0,0, (0)1d xdx

26、dxxxxdtdtdt121222121(0)1,(0)0(1)dyyydtdyyyyydt %解算微分方程解算微分方程tspan=0,30; y0=1;0;tt,yy=ode45(DyDt,tspan,y0); figure(1)plot(tt,yy(:,1)xlabel(t),title(x(t)dtdxyxy21,解解:令令，上式寫成一階微分方程組方式，上式寫成一階微分方程組方式%畫相平面圖畫相平面圖(函數(shù)和其導(dǎo)函數(shù)和其導(dǎo)數(shù)勾畫的曲線稱為數(shù)勾畫的曲線稱為相軌相軌跡跡)figure(2)plot(yy(:,1),yy(:,2)xlabel(位移位移),ylabel(速度速度) functi

27、on ydot=DyDt(t,y)mu=2; ydot=y(2);mu*(1-y(1)2)*y(2)-y(1);據(jù)以上方程組，編寫據(jù)以上方程組，編寫M函數(shù)文件函數(shù)文件DyDt.m051015202530-2.5-2-1.5-1-0.500.511.522.5tx(t)-2.5-2-1.5-1-0.500.511.522.5-4-3-2-101234位 移速度4.2矩陣和代數(shù)方程矩陣和代數(shù)方程4.2.1 矩陣運(yùn)算和特征參數(shù)矩陣運(yùn)算和特征參數(shù)4.2.2 矩陣的變換和特征值分解矩陣的變換和特征值分解4.2.3 線性方程的解線性方程的解4.2.4 普通代數(shù)方程的解普通代數(shù)方程的解4.2.1 矩陣運(yùn)算和

28、特征參數(shù)矩陣運(yùn)算和特征參數(shù) 矩陣與標(biāo)量之間的四那么運(yùn)算與數(shù)組運(yùn)算一樣矩陣與標(biāo)量之間的四那么運(yùn)算與數(shù)組運(yùn)算一樣 矩陣和矩陣之間的四那么運(yùn)算矩陣和矩陣之間的四那么運(yùn)算 矩陣和矩陣之間的加減運(yùn)算與數(shù)組運(yùn)算一樣矩陣和矩陣之間的加減運(yùn)算與數(shù)組運(yùn)算一樣 設(shè)設(shè) A 是一個(gè)是一個(gè) mn 矩陣，矩陣，B 是一個(gè)是一個(gè) pq 矩陣，當(dāng)矩陣，當(dāng) np 時(shí)，兩個(gè)矩陣可以相乘，乘積為時(shí)，兩個(gè)矩陣可以相乘，乘積為 mq 矩陣。矩陣乘矩陣。矩陣乘法不可逆。在法不可逆。在 MATLAB 中，矩陣乘法由中，矩陣乘法由“*實(shí)現(xiàn)。實(shí)現(xiàn)。 矩陣除法在實(shí)踐中主要用于求解線性方程組矩陣除法在實(shí)踐中主要用于求解線性方程組 矩陣轉(zhuǎn)置：矩陣轉(zhuǎn)

29、置： 符號符號“ 實(shí)現(xiàn)矩陣的轉(zhuǎn)置操作。對于實(shí)數(shù)矩陣，實(shí)現(xiàn)矩陣的轉(zhuǎn)置操作。對于實(shí)數(shù)矩陣， “ 表示矩陣轉(zhuǎn)置，對于復(fù)數(shù)矩陣，表示矩陣轉(zhuǎn)置，對于復(fù)數(shù)矩陣，“ 實(shí)現(xiàn)共軛轉(zhuǎn)實(shí)現(xiàn)共軛轉(zhuǎn)置。對于復(fù)數(shù)矩陣，假設(shè)想要實(shí)現(xiàn)非共軛轉(zhuǎn)置，可以置。對于復(fù)數(shù)矩陣，假設(shè)想要實(shí)現(xiàn)非共軛轉(zhuǎn)置，可以運(yùn)用符號運(yùn)用符號“ . 。format rat %有理格式顯示有理格式顯示A=magic(2) + j*pascal(2) A = 1 + 1i 3 + 1i 4 + 1i 2 + 2iA1=A A1=1 - 1i 4 - 1i %共軛轉(zhuǎn)置共軛轉(zhuǎn)置 3 - 1i 2 - 2iA2=A. A2=1 + 1i 4 + 1i %非共軛轉(zhuǎn)置，

30、數(shù)組運(yùn)算操非共軛轉(zhuǎn)置，數(shù)組運(yùn)算操作作 3 + 1i 2 + 2i例例4.2-2 矩陣和數(shù)組轉(zhuǎn)置操作的差別矩陣和數(shù)組轉(zhuǎn)置操作的差別計(jì)算矩陣標(biāo)量特征參數(shù)計(jì)算矩陣標(biāo)量特征參數(shù)-秩秩,跡跡,行列式的指令行列式的指令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=t

31、race(A)【例【例4.2-3】矩陣標(biāo)量特征參數(shù)計(jì)算例如?！烤仃嚇?biāo)量特征參數(shù)計(jì)算例如。A = 1 4 7 2 5 8 3 6 9r = 2d3 =0d2 = -3t = 15 【例【例4.2-4】矩陣標(biāo)量特征參數(shù)的性質(zhì)。】矩陣標(biāo)量特征參數(shù)的性質(zhì)。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) %任何符合矩陣乘法規(guī)那么的兩個(gè)矩陣乘任何符合矩陣乘法規(guī)那么的兩個(gè)矩陣乘積積tBA=trace(B*A) %的的“跡不變。同階乘積跡不變。同階乘積“跡不變跡不變tC

32、D=trace(C*D) %兩個(gè)兩個(gè)“內(nèi)維相等矩陣的乘積內(nèi)維相等矩陣的乘積“跡不變跡不變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 非同階矩陣乘積行列非同階矩

33、陣乘積行列式不等于各矩陣行列式不等于各矩陣行列式之乘積式之乘積udx=diff(X) %計(jì)算向量計(jì)算向量X的前向差分的前向差分uFX=gradient(F) %求一元求一元(函數(shù)函數(shù))梯度梯度uFX, FY =gradient(F) %求二元函數(shù)梯度求二元函數(shù)梯度4.1.1 近似數(shù)值極限及導(dǎo)數(shù)近似數(shù)值極限及導(dǎo)數(shù)S1=sum(X,1) % sum按列向求和得按列向求和得(1n)數(shù)組數(shù)組, =sum(X) (X多行多行)S2=sum(X,2) % sum按行向求和得按行向求和得(n1)數(shù)組數(shù)組Scs=cumsum(X)%沿沿X列向求累計(jì)和列向求累計(jì)和, 仍是數(shù)組仍是數(shù)組, 最后一行最后一行等于等于

34、Sx4.1.2 數(shù)值求和與近似數(shù)值積分?jǐn)?shù)值求和與近似數(shù)值積分St=trapz(t,X) 或或 St=dt*trapz(X) %梯形法求積梯形法求積分分Sct=cumtrapz(t,X) 或或 Sct=dt*cumtrapz(X) %梯形法沿列向求梯形法沿列向求X關(guān)于關(guān)于x的累計(jì)積分的累計(jì)積分,最后一最后一個(gè)值等于個(gè)值等于StX= 1 2; 3 4 S1= 4 6 S2= 3 7sum(S1) = 10sum(S1,1) = 4 6Review4.1.3 計(jì)算精度可控的數(shù)值積分計(jì)算精度可控的數(shù)值積分一重積分一重積分(quadrature精度可控精度可控)：S1=quad(fun,a,b,tol)

35、 %自順應(yīng)自順應(yīng)Simpson法法S2=quadl(fun,a,b,tol) %自順應(yīng)羅巴托自順應(yīng)羅巴托 Lobatlo法法二重積分二重積分(精度可控精度可控) ：S3=dblquad(fun,xmin,xmax,ymin,ymax,tol)三重積分三重積分(精度可控精度可控) ： S4=triplequad(fun,xmin,xmax,ymin,ymax, zmix,zmax,tol) fun:可為字符串可為字符串,內(nèi)聯(lián)對象內(nèi)聯(lián)對象,匿名函數(shù)匿名函數(shù),函數(shù)句柄函數(shù)句柄,留意乘除留意乘除法運(yùn)算一定加點(diǎn)法運(yùn)算一定加點(diǎn).用數(shù)組運(yùn)算用數(shù)組運(yùn)算 tol: 默許積分的絕對精度為默許積分的絕對精度為10-

36、6Review4.1.4 函數(shù)極值的數(shù)值求解函數(shù)極值的數(shù)值求解 x,fval,exitflag,output=fminbnd(fun,x1,x2,options) %一元函數(shù)在區(qū)間一元函數(shù)在區(qū)間bound(x1,x2)中極小值中極小值 x,fval: 極值點(diǎn)坐標(biāo)和對應(yīng)目的函數(shù)極值極值點(diǎn)坐標(biāo)和對應(yīng)目的函數(shù)極值x,fval,exitflag,output=fminsearch(fun,x0,options) %單純形法求搜索起點(diǎn)單純形法求搜索起點(diǎn)x0附近多元函數(shù)極值點(diǎn)附近多元函數(shù)極值點(diǎn) % x每列代表一個(gè)候選極值點(diǎn)，各列按目的函數(shù)每列代表一個(gè)候選極值點(diǎn)，各列按目的函數(shù)極小值遞增順序極小值遞增順序

37、%x(:,1)對應(yīng)的目的函數(shù)極小值點(diǎn)由對應(yīng)的目的函數(shù)極小值點(diǎn)由fval給出給出 fun: 字符串字符串,內(nèi)聯(lián)函數(shù)內(nèi)聯(lián)函數(shù),匿名函數(shù)匿名函數(shù),函數(shù)句柄函數(shù)句柄, 留意留意乘除法運(yùn)算一定加點(diǎn)乘除法運(yùn)算一定加點(diǎn).用數(shù)組運(yùn)算用數(shù)組運(yùn)算 options: 配置優(yōu)化參數(shù)配置優(yōu)化參數(shù),可略可略 exitflag: 給出大于給出大于0的數(shù)的數(shù),那么勝利搜索到極值點(diǎn)那么勝利搜索到極值點(diǎn) output: 給出詳細(xì)的優(yōu)化算法和迭代次數(shù)給出詳細(xì)的優(yōu)化算法和迭代次數(shù)Reviewrank(A) %求秩求秩Rankdet(A) %求行列式求行列式Determinanttrace(A) %求跡求跡(Trace)，即矩陣主對角

38、元素的和，即矩陣主對角元素的和Reviewt,Y=ode45(odefun,tspan,y0) % 4階龍格庫塔階龍格庫塔 數(shù)數(shù)值法值法odefun: 待求解一階微分方程組的函數(shù)文件句柄待求解一階微分方程組的函數(shù)文件句柄tspan: 自變量微分二元區(qū)間自變量微分二元區(qū)間t0, tfy0: 一階微分方程組的一階微分方程組的(n*1)初值列向量初值列向量ode45運(yùn)用方法運(yùn)用方法:t-二元區(qū)間的點(diǎn)系列二元區(qū)間的點(diǎn)系列Y-原函數(shù)在微分區(qū)間點(diǎn)系列上的函數(shù)值原函數(shù)在微分區(qū)間點(diǎn)系列上的函數(shù)值4.1.5 常微分方程常微分方程Ordinary Differential Equation的的數(shù)值解數(shù)值解4.2.

39、1 矩陣運(yùn)算和特征參數(shù)矩陣運(yùn)算和特征參數(shù)4.2.2 矩陣的變換和特征值分解矩陣的變換和特征值分解R, ci=rref(A) %借助初等變換將借助初等變換將A縮減行變成行階梯縮減行變成行階梯矩陣矩陣R。B=orth(A)可得到矩陣可得到矩陣A的正交基，的正交基，B的列與的列與A的列可張成一樣的列可張成一樣的空間，而且的空間，而且B的列是正交的，因此的列是正交的，因此B*B=eye(rank(A)，B的的列數(shù)正好是列數(shù)正好是A的秩。的秩。 X=null(A) %A矩陣零空間的全部正交基，滿足矩陣零空間的全部正交基，滿足AX=0, X*X=I。B=orth(A)V, D=eig(A) %A矩陣的特征

40、值、特征向量分解，使矩陣的特征值、特征向量分解，使AV=VD。ci 是行數(shù)組是行數(shù)組, 其元素表示其元素表示A中線性獨(dú)立列的序號。中線性獨(dú)立列的序號。length(ci)=rank(A)【例【例4.2-5】行階梯陣簡化指令】行階梯陣簡化指令rref計(jì)算結(jié)果的含義。計(jì)算結(jié)果的含義。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

41、(:,1:3)*R(1:3,4) %A前三列線形組前三列線形組合合err=norm(A(:,4)-aa) aa = 13 8 12 1err = 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 矩陣零空間及含義。設(shè) 是矩陣 的零空間,即n lXm nA0,m nn lm lAX()m nrank

42、 Anl有n=size(A,2)l=size(X,2)Rank_A=rank(A)n-l=rank(A)【例【例4.2-7】簡單實(shí)陣的特征分解?！亢唵螌?shí)陣的特征分解。eig，cdf2rdf 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=

43、cdf2rdf(V,D) VR =0.7746 0 0.0430 -0.6310DR =0.8333 2.4438 -2.4438 0.8333 %3A1=V*D/V %因計(jì)算誤差有很因計(jì)算誤差有很小虛部小虛部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.

44、6667err1 =6.7532e-016err2 =4.4409e-016 Frobenius矩陣范數(shù)矩陣范數(shù) 實(shí)數(shù)塊實(shí)數(shù)塊對角陣對角陣實(shí)陣實(shí)陣矩陣矩陣Frobenius范數(shù)范數(shù),類似向量類似向量2范數(shù)，不同于范數(shù)，不同于矩陣矩陣2范數(shù)范數(shù)4.2.2 線性方程的解線性方程的解對于方程組對于方程組Amnx=b (m個(gè)方程，個(gè)方程，n個(gè)未知數(shù)個(gè)未知數(shù))， 當(dāng)向量當(dāng)向量b在矩陣在矩陣A列向量所張空間中列向量所張空間中, rank(A,b) =rank(A)= r a) 假設(shè)假設(shè)n =r，那么解獨(dú)一。，那么解獨(dú)一。 b) 假設(shè)假設(shè)n r，那么解不獨(dú)一。，那么解不獨(dú)一。 當(dāng)向量當(dāng)向量b不在矩陣不在矩陣

45、A列向量所張空間中列向量所張空間中, 那么無準(zhǔn)確解那么無準(zhǔn)確解但存在最小二乘解。但存在最小二乘解。1. 線性方程解的普通討論線性方程解的普通討論 matlab定義的左除運(yùn)算可以很方便地解上述方程組：定義的左除運(yùn)算可以很方便地解上述方程組： x=Ab (x=A-1b只適用于只適用于A非奇特時(shí)非奇特時(shí))2.除法運(yùn)算解方程解除法運(yùn)算解方程解 當(dāng)當(dāng)m=n時(shí)時(shí), “恰定方程恰定方程 當(dāng)當(dāng)mn時(shí)時(shí), “超定方程超定方程 當(dāng)當(dāng)mr, 解不獨(dú)一解不獨(dú)一% 3求特解和通解，并對由之構(gòu)成的全解進(jìn)展驗(yàn)算求特解和通解，并對由之構(gòu)成的全解進(jìn)展驗(yàn)算xs=Ab; xg=null(A); % xg是齊次方程是齊次方程Ax=0

46、的解的解c=rand(1); ba=A*(xs+c*xg) % ba是是A乘乘 “一個(gè)隨機(jī)的全解一個(gè)隨機(jī)的全解Warning: Rank deficient, rank = 2, tol = 1.8757e-014.ba =13.0000 14.0000 15.0000 16.0000 例例4.2-8: 求方程求方程 的解的解15913261014371115481216xnorm(ba-b) ans = 1.4e-014 例例4.2-9: “逆陣法和左除法解恰定方程的性能對逆陣法和左除法解恰定方程的性能對比。比。測試陣測試陣產(chǎn)生指定異常值和特殊帶寬的隨機(jī)陣產(chǎn)生指定異常值和特殊帶寬的隨機(jī)陣3.

47、 矩陣的逆矩陣的逆A_1=inv(A) %求非奇特陣求非奇特陣A的逆的逆randn(state,0); A=gallery(randsvd,300,2e13); %產(chǎn)生條件數(shù)產(chǎn)生條件數(shù) 為為2e13的的300階隨機(jī)矩陣階隨機(jī)矩陣;x=ones(300,1); %定義真解定義真解b=A*x;cond(A) %驗(yàn)算矩陣條件數(shù)驗(yàn)算矩陣條件數(shù), 結(jié)果結(jié)果a1.9978e+013, 值越大值越大, 陣越病態(tài)陣越病態(tài),%求逆法求逆法tic; xi=inv(A)*b; ti=toceri=norm(x-xi) rei=norm(A*xi-b)/norm(b) %左除法左除法tic; xd=Ab;td=toc

48、erd=norm(x-xd)red=norm(A*xd-b)/norm(b)ti =0.0185eri =0.0883rei =0.0051td =0.0066erd =0.0298red =8.7810e-015矩陣計(jì)算對于誤差越敏感矩陣計(jì)算對于誤差越敏感, 數(shù)值穩(wěn)定性差數(shù)值穩(wěn)定性差 求解恣意函數(shù)求解恣意函數(shù)f(x)=0(能夠無解能夠無解,單解單解,多解多解)的步驟的步驟:作圖獲取初步近似解作圖獲取初步近似解 察看察看f(x)與橫軸的交點(diǎn)坐標(biāo)，用與橫軸的交點(diǎn)坐標(biāo)，用zoom放大，用放大，用ginput得較準(zhǔn)確些的交點(diǎn)坐標(biāo)值。得較準(zhǔn)確些的交點(diǎn)坐標(biāo)值。4.3 普通代數(shù)方程的解普通代數(shù)方程的解 用

49、泛函指令求準(zhǔn)確解用泛函指令求準(zhǔn)確解 x,favl = fzero(fun,x0) %一元函數(shù)求零點(diǎn)一元函數(shù)求零點(diǎn) x,fval = fsolve(fun,x0) %解非線性方程組解非線性方程組 fun: 內(nèi)聯(lián)對象內(nèi)聯(lián)對象,匿名函數(shù)匿名函數(shù),函數(shù)句柄函數(shù)句柄,字符串字符串;被解函被解函數(shù)自變量普通用數(shù)自變量普通用x, 留意乘除法運(yùn)算一定加點(diǎn)留意乘除法運(yùn)算一定加點(diǎn).用數(shù)組用數(shù)組運(yùn)算運(yùn)算 x0: 零點(diǎn)初始猜測值零點(diǎn)初始猜測值. x0為標(biāo)量時(shí)取與之最接近的為標(biāo)量時(shí)取與之最接近的零點(diǎn)零點(diǎn); x0取取a,b時(shí)時(shí),在此區(qū)間內(nèi)尋覓一個(gè)零點(diǎn)在此區(qū)間內(nèi)尋覓一個(gè)零點(diǎn). x: 所求零點(diǎn)的自變量值所求零點(diǎn)的自變量值 f

50、val: 函數(shù)值函數(shù)值-2-1.5-1-0.500.511.522.5-0.6-0.4-0.200.20.40.6ty(t)例4.2-10: 求f(t)=(sin2t)e-0.1t0.5|t|的零點(diǎn)-10-50510-5-4-3-2-101ty(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.

51、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 概率分布和統(tǒng)計(jì)分析概率分布和統(tǒng)計(jì)分析4.3.1 概率函數(shù)、分布函數(shù)、逆分布概率函數(shù)、分布函數(shù)、逆分布函數(shù)和隨機(jī)數(shù)的發(fā)生函數(shù)和隨機(jī)數(shù)的發(fā)生4.3.2 隨機(jī)數(shù)發(fā)生器和隨機(jī)數(shù)發(fā)生器和 統(tǒng)計(jì)分析指統(tǒng)計(jì)分析指令令4.3 概率分布和統(tǒng)計(jì)分析概率分布和統(tǒng)計(jì)分析1. 二項(xiàng)分布二項(xiàng)分布Binomial distrib

52、ution 4.3.1 概率函數(shù)、分布函數(shù)、逆分布函數(shù)和隨機(jī)數(shù)的發(fā)生pk=binopdf(k,N,p) 事件事件A發(fā)生發(fā)生k次的概率次的概率(1)kkn knP XkC ppFk=binocdf(k,N,p) 事件事件A發(fā)生次數(shù)不大于發(fā)生次數(shù)不大于k的的概率概率0( )(1)kjjnjnjF kC ppR=binornd(N,p) 產(chǎn)生符合二項(xiàng)分布B(N,p)的(mn) 隨機(jī)數(shù)組元素值為事件A能夠發(fā)生的次數(shù)。Binomial Cumulative Distribution FunctionBinomial Probability Density FunctionN=100;p=0.5;k=0:

53、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),Children),Color,r,Marker,+,MarkerSize,4)set(get(h(2),Ylabel),String,cdf)xlabel(k) grid on 010203040506070809010000.050.1010203040506070809010000.51例例4.3-1. 畫出畫出N=100，p=0.5情況下的二項(xiàng)分布概率特性情況下的二項(xiàng)分布概率特性曲線。曲線。4.3.1 概率函數(shù)、分布函數(shù)、逆分布函數(shù)和隨機(jī)數(shù)的發(fā)生010203040506070809010000.050.1kpdf010203040506070809010000.51cdfpx=nor