ch4數(shù)值計(jì)算2

1、第 4 章 數(shù)值計(jì)算本章主要講的函數(shù)有1、 近似數(shù)值極限及導(dǎo)數(shù)的函數(shù)有diffgradient2 、 數(shù)值求和與近似數(shù)值積分函數(shù)有sumsumsumtrapzcumtranpz3、 計(jì)算精度可控的數(shù)值積分integral4、 函數(shù)極值的數(shù)值求解fminbndfminserch5、 常微分方程的數(shù)值求解ode456、 矩陣和代數(shù)方程求矩陣的秩：rank跡：trace行列式：det階梯矩陣變換：rrefA矩陣零空間的全部正交基：nullA矩陣值空間的全部正交基：orthA矩陣的特征值，特征向量分解：eig7、 一般代數(shù)方程的解一元函數(shù)零點(diǎn)指令：fzero解非線性方程組：fsolve4.1 數(shù)值微積

2、分（1）重視有限精度浮點(diǎn)表示的離散本質(zhì)，不要貿(mào)然自行編寫數(shù)值計(jì)算程序進(jìn)行求導(dǎo)和求極限運(yùn)算（2）數(shù)值導(dǎo)數(shù)一定要求取，自變量的增量大于原數(shù)據(jù)精度的10倍以上。（3）解極限值，積分，微分方程數(shù)值計(jì)算時(shí)，要盡量使用MATLAB提供的指令，嚴(yán)格遵守指令的使用規(guī)則。特別說明：MATLAB沒有提供數(shù)值求導(dǎo)和極限的指令。其中概念法求差分dx=diff(X)求差分FX=gradient(F)FX,FY=gradient(F)4.1.1 近似數(shù)值極限及導(dǎo)數(shù)在數(shù)值法近似求極限和理論有時(shí)可能不一致（舉例說明）,數(shù)值法求極限輕易不用【例4.1-1】設(shè)試用機(jī)器零閾值eps替代理論0計(jì)算極限理論分析：%不可信的極限的數(shù)值近

3、似計(jì)算x=eps;L1=(1-cos(2*x)/(x*sin(x),%L2=sin(x)/x,% L1 = 0L2 = 1 %可信的極限的符號(hào)計(jì)算syms 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 【例4.1-2】已知x=sin(t)，求該函數(shù)在區(qū)間0 2pi中的近似導(dǎo)函數(shù)。本例考察的是，在將連續(xù)信號(hào)進(jìn)行離散化時(shí)，采樣間隔影響導(dǎo)函數(shù)。（1）計(jì)算數(shù)值導(dǎo)數(shù)時(shí)，自變量的額增量取得過?。ㄔ趀ps數(shù)量級(jí)）d=pi/100;t=0:d:2*pi;x=sin(t);dt=5*ep

4、s;%增量為epsx_eps=sin(t+dt);dxdt_eps=(x_eps-x)/dt;%以dt=5*eps為增量算得的數(shù)值導(dǎo)數(shù)plot(t,x,'LineWidth',5)hold on plot(t,dxdt_eps)hold offlegend('x(t)','dx/dt')xlabel('t') 圖 4.1-1 增量過小引起有效數(shù)字嚴(yán)重丟失后的毛刺曲線（2）計(jì)算數(shù)值導(dǎo)數(shù)時(shí)，自變量的增量取得適當(dāng)x_d=sin(t+d);dxdt_d=(x_d-x)/d;%以d=pi/100為增量得的數(shù)值導(dǎo)數(shù)plot(t,x,'

5、LineWidth',5)hold onplot(t,dxdt_d)hold offlegend('x(t)','dx/dt')xlabel('t') 圖 4.1-2 增量適當(dāng)所得導(dǎo)函數(shù)比較光滑現(xiàn)象差距合理解釋如果dt步長(zhǎng)太小，f(t+dt)與f(t)數(shù)值十分接近，它們的高維有效數(shù)字完全相同。f(t+dt)與f(t)差高位有效數(shù)字消失，導(dǎo)致精度急劇下降?！纠?.1-3】已知x=sin(t),采用diff和gradient計(jì)算該函數(shù)在區(qū)間0 2*pi中的近似導(dǎo)數(shù)clfd=pi/100;%t=0:d:2*pi;x=sin(t);dxdt_di

6、ff=diff(x)/d;%dxdt_grad=gradient(x)/d;%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',

7、'Location','North')xlabel('t'),box offhold offsubplot(1,2,2)kk=(length(t)-10):length(t);%thold onplot(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'

8、;,'dxdt_diff','Location','SouthEast')xlabel('t'),box offhold off 圖 4.1-3 diff和gradient求數(shù)值近似導(dǎo)數(shù)的異同比較4.1.2 數(shù)值求和與近似數(shù)值積分Sx=sum(X);沿著列方向求和Scs=cumsum(X)：沿著列方向求累計(jì)和St=trapz(x,y)采用梯形法沿著類方向求函數(shù)y關(guān)于自變量x的積分Sct=cumtrapz(x,y):采用梯形沿著列方向求函數(shù)y關(guān)于自變量x的累計(jì)積分假如Xm*n數(shù)組，sum（X）的計(jì)算結(jié)果表示第K列全體元素之和【例 4

9、.1-4】 求積分其中y=0.2+sintcleard=pi/8;%分區(qū)間的區(qū)間間隔t=0:d:pi/2;%包含5個(gè)采樣點(diǎn)的一維數(shù)組y=0.2+sin(t);%5個(gè)點(diǎn)處的函數(shù)數(shù)值數(shù)組s=sum(y);%求出所有函數(shù)采樣值之和s_sa=d*s;%高度為函數(shù)采樣值得素有小矩形面積之和 <6>s_ta=d*trapz(y);%連接個(gè)函數(shù)采樣值得折線下的面積<7>disp('sum求得積分',blanks(3),'trapz求得積分')disp(s_sa, s_ta)t2=t,t(end)+d;%因采用stairs繪圖需要而寫y2=y,nan;%

10、因采用stairs繪圖需要而寫stairs(t2,y2,':k')%用虛線下面積表示d*sum的幾何意義hold onplot(t,y,'r','LineWidth',3)%用紅色折線下面積表示d*trapz幾何意義h=stem(t,y,'LineWidth',2);%用藍(lán)空心桿線表示函數(shù)采樣值set(h(1),'MarkerSize',10)axis(0,pi/2+d,0,1.5)%使橫坐標(biāo)恰好是0,5*dhold offshg sum求得積分 trapz求得積分 1.5762 1.3013圖 4.1-4 sum

11、 和trapz求積模式示意4.1.3 計(jì)算精度可控的數(shù)值積分4.1.2數(shù)值積分中是根據(jù)傳統(tǒng)的方法編寫了2個(gè)指令，用以計(jì)算數(shù)值積分。但由于實(shí)際工程中以上2個(gè)指令可能存在精度難以控制，或不能處理廣義積分，還可能運(yùn)算速度慢。所以根據(jù)最新算法設(shè)計(jì)了2個(gè)函數(shù)。q=integral(fun,xmin,xmax)意義：在xmin到xmax區(qū)間上計(jì)算fun函數(shù)的精確積分（簡(jiǎn)單調(diào)用）q=integral(fun,xmin,xmax,Name,Value)意義：在xmin到xmax區(qū)間上計(jì)算fun函數(shù)的精確積分（復(fù)雜調(diào)用）說明：fun:表示被積函數(shù)的M碼匿名函數(shù)或函數(shù)句柄 【例 4.1-5】（1）x=linspa

12、ce(0.01,1.2,50);構(gòu)造50個(gè)自變量X數(shù)組g1=x.(-0.2);g2=x.5;計(jì)算兩個(gè)函數(shù)對(duì)應(yīng)自變量的函數(shù)值plot(x,g1,'-r.',x,g2,'-b*')axis(0,1.2,0,3)legend('g_1(x)=1/x0.2','g_2(x)=x5','Location','North')title('x位于0,1間的g_1(x)曲線與g_2(x)曲線所夾的區(qū)域') 圖 4.1-5 待求面積的區(qū)域（2）采用標(biāo)量型匿名函數(shù)計(jì)算積分format longG1=(x

13、)x.-0.2;%<8>構(gòu)造g1(x)匿名函數(shù)Q1=integral(G1,0,1)%<9>在0 1區(qū)間計(jì)算g1(x)曲線下的面積G2=(x)x.5;%<10>構(gòu)造g2(x)匿名函數(shù)Q2=integral(G2,0,1)%<11>在0 1區(qū)間計(jì)算g2(x)曲線下的面積S12=Q1-Q2 %<12> 兩曲線所夾區(qū)域的面積Q1 = 1.250000027856048Q2 = 0.166666666666667S12 = 1.083333361189381 （3）采用陣列型匿名函數(shù)計(jì)算積分G=(x)x.-0.2;5;%<13>構(gòu)

14、造陣列性匿名函數(shù)Q=integral(G,0,1,'ArrayValued',true)%<14>復(fù)雜調(diào)用格式的 積分結(jié)果S=1,-1*Q% <15> Q = 1.250000027856048 0.166666666666667S = 1.083333361189381 （4）符號(hào)積分驗(yàn)證syms t%定義符號(hào)變量Gsym=vpa(int(t.-0.2;5,0,1);%計(jì)算兩個(gè)函數(shù)具有32精度的積分值Ssym=Gsym(1)-Gsym(2)%<17>至少31位精度的曲線所圍區(qū)域面積Ssym =1.0833333333333333333333

15、333333333 4.1.4 函數(shù)極值的數(shù)值求解只有處理極小值（局域極值）指令，沒有專門的極大值指令（極大值是極小值f(x)的負(fù)數(shù)即-f(x)）。x,fval,exitflag,output=fminbnd(fun,x1,x2,options)求一元函數(shù)在區(qū)間（x1,x2）中極小值x,fval,exitflag,output=fminsearch(fun,x0,options)單純形法求多遠(yuǎn)函數(shù)極值點(diǎn)說明：(1) fminbnd: 輸入?yún)?shù)說明 fun:兩個(gè)函數(shù)中都有的fun:待求目標(biāo)函數(shù)書寫方法 M文件的函數(shù)句柄 字符串 內(nèi)聯(lián)對(duì)象 匿名函數(shù)x1左邊界；x2：右邊界options:配置優(yōu)化參

16、數(shù)（通常默認(rèn)即不用書寫）輸出參數(shù)說明x,極值點(diǎn)fval,目標(biāo)函數(shù)exitflag,若大于0的數(shù)，說明成功的搜索島極值點(diǎn)output：給出優(yōu)化算法和迭代次數(shù) （2）fminsearch輸入?yún)?shù)說明fun：同fminbndx0,搜索七點(diǎn)的向量或一組搜索起點(diǎn)的矩陣 當(dāng)采用單個(gè)起點(diǎn)搜索時(shí)，輸出量x也是一個(gè)單點(diǎn) 當(dāng)采用多個(gè)搜索起點(diǎn)（矩陣）時(shí)，該矩陣的每一列代表一個(gè)候選極值點(diǎn)。 搜索到的極值點(diǎn)按照目標(biāo)函數(shù)遞增次序排列。極值點(diǎn)x(:,1)對(duì)應(yīng)的目標(biāo)函數(shù)極小值由f val給出。fval,目標(biāo)函數(shù)（同fminbnd）exitflag,若大于0的數(shù)，說明成功的搜索到極值點(diǎn)（同fminbnd）輸出參數(shù)說明x,極值點(diǎn)

17、（同fminbnd）fval,目標(biāo)函數(shù)（同fminbnd）exitflag,若大于0的數(shù)，說明成功的搜索島極值點(diǎn)（同fminbnd）output：給出優(yōu)化算法和迭代次數(shù) （同fminbnd）【例4.1-7】已知區(qū)間，求函數(shù)的最小值。本題目的思路：（1）符號(hào)計(jì)算求極值的局限性（2）fminbnd求極小值的局限性（3）求最小值時(shí)，需要整個(gè)區(qū)間的函數(shù)信息和圖形法功用（先畫整個(gè)區(qū)間的圖示，限定好區(qū)間，再求，不能輕信fminbnd函數(shù)的結(jié)果）（1）用“導(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ù)xs

18、0=solve(yd,x)%求導(dǎo)函數(shù)為0的自變量值xs0<4>yd_xs0=vpa(subs(yd,x,xs0),6)%驗(yàn)算用：導(dǎo)函數(shù)在xs0處為0嗎y_xs0=vpa(subs(y,x,xs0),6)%在xs0處的函數(shù)值xs0 =0.050838341410271656880659496266968yd_xs0 =2.29589*10(-41)y_xs0 =-0.00126332 （2）采用優(yōu)化算法求極小值x1=-10;x2=10;%搜索區(qū)間的邊界yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);%采用匿名函數(shù)形式定義被求極小值的函數(shù)y(

19、x)xn0,fval,exitflag,output=fminbnd(yx,x1,x2)%<9>% xn0,fva1分別是極值點(diǎn)和函數(shù)極小值xn0 = 2.514797840754235fval = -0.499312445280039exitflag = 1output = iterations: 13 funcCount: 14 algorithm: 'golden section search, parabolic interpolation' message: 1x112 char （3）繪圖觀察最小值xx=-10:pi/200:10;%采用點(diǎn)應(yīng)做夠密yxx

20、=subs(y,x,xx);plot(xx,yxx)xlabel('x'),grid on 圖 4.1-5 在-10, 10區(qū)間中的函數(shù)曲線（4）距圖形觀察，重設(shè)fminbnd的搜索區(qū)間x11=6;x2=10;%搜索區(qū)間的邊界yx=(x)(sin(x)2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);%采用匿名函數(shù)形式定義被求極小值的函數(shù)y(x)xn00,fval,exitflag,output=fminbnd(yx,x11,x2)%<16>% xn0,fva1分別是極小值點(diǎn)和函數(shù)極小值xn00 = 8.023562824723015fval = -

21、3.568014059128578exitflag = 1output = iterations: 9 funcCount: 10 algorithm: 'golden section search, parabolic interpolation' message: 1x112 char 【例4.1-8】求的極小值點(diǎn)。它既是注明的Rosenbrock's 'Banna'測(cè)試函數(shù)，它的理論極小值是x=1,y=1（1）本例采用匿名函數(shù)表示測(cè)試函數(shù)（含有y,x變量的自變量，寫成x(1),x(2)）ff=(x)(100*(x(2)-x(1).2)2+(1-x(

22、1)2); （2）采用單純形法求極小值（sx,一組極小值點(diǎn)坐標(biāo)，sfval目標(biāo)值）format short gx0=-5,-2,2,5;-5,-2,2,5;%提供4個(gè)搜索起點(diǎn)sx,sfval,sexit,soutput=fminsearch(ff,x0)%sx給出一組使優(yōu)化函數(shù)值非減的局部極小點(diǎn) sx = 0.99998 -0.68971 0.41507 8.0886 0.99997 -1.9168 4.9643 7.8004sfval = 2.4112e-010sexit = 1soutput = iterations: 384 funcCount: 615 algorithm: '

23、Nelder-Mead simplex direct search' message: 1x196 char （3）檢查目標(biāo)函數(shù)format short edisp(ff(sx(:,1),ff(sx(:,2),ff(sx(:,3),ff(sx(:,4) 2.4112e-010 5.7525e+002 2.2967e+003 3.3211e+005 4.1.5 常微分方程的數(shù)值解MATLAB提供一套好的求微分方程的指令t,Y=ode45(odefun,tspan,y0)采用4階Runge-Kutta數(shù)值積分法解微分方程說明：（1） 輸入變量：odefun:是待解微分方程的函數(shù)文件句柄。要

24、求函數(shù)文件的輸出形式必須是待解函數(shù)的一階導(dǎo)數(shù)。也就是不管原方程是不是一階導(dǎo)數(shù)，必須化成一階導(dǎo)數(shù)的形式才可以用ode45求解。一解微分方程形式為：,y是（n*1）向量tspan：常被賦成二元向量t0,tf,此時(shí)tspan用來(lái)定義求數(shù)值解的時(shí)間區(qū)間。y0： 一解微分方程組的n*1初值列向量t： 所求數(shù)值解的自變量數(shù)據(jù)列向量（假定其數(shù)據(jù)長(zhǎng)度為N），而Y則是（N*n）矩陣，輸出量Y行中第K列Y(:,K)，就是微分方程y的第k分量的解【例4.1-9】求解微分方程，在初始條件，情況下的解，并圖示（1） 把高階微分方程改寫成一階微分方程組，令寫成微分方程系數(shù)矩陣（雅可比系數(shù)陣）的形式，=2（2） 根據(jù)上述一

25、階微分方程組編寫M函數(shù)文件DyDt.mfunction ydot=DyDt(t,y)mu=2ydot=y(2);mu*(1-y(1)2)*y(2)-y(1);注意一階導(dǎo)數(shù)ydot是2*1列向量（3） 解微分方程tspan=0,30;%求解的時(shí)間區(qū)間y0=1;0;%初值向量應(yīng)與DyDt.m文件中y形式一致tt,yy=ode45(DyDt,tspan,y0);%<3>plot(tt,yy(:,1)xlabel('t'),title('x(t)') 圖 4.1-6 微分方程解（4）plot(yy(:,1),yy(:,2)%函數(shù)和其導(dǎo)函數(shù)勾畫的曲線稱為“相軌

26、跡”xlabel('位移'),ylabel('速度') 圖 4.1-7 平面相軌跡4.2 矩陣和代數(shù)方程4.2.1 矩陣的標(biāo)量特征參數(shù)【例 4.2-1】A=reshape(1:9,3,3)%r=rank(A)%d3=det(A)%d2=det(A(1:2,1:2)%t=trace(A)% A = 1 4 7 2 5 8 3 6 9 r = 2 d3 = 0 d2 = -3 t = 15 【例4.2-2】format short%rng default%A=rand(3,3);%B=rand(3,3);%C=rand(3,4);D=rand(4,3);tAB=tr

27、ace(A*B)%tBA=trace(B*A) tCD=trace(C*D)%tDC=trace(D*C) tAB = 3.7479tBA = 3.7479tCD = 3.3399tDC = 3.3399 d_A_B=det(A)*det(B)dAB=det(A*B)dBA=det(B*A) d_A_B = -0.0852dAB = -0.0852dBA = -0.0852 dCD=det(C*D)dDC=det(D*C)% dCD = -0.0557dDC = 1.5085e-017 4.2.2 矩陣的變換和特征值分解【例4.2-3】（1）A=magic(4)%R,ci=rref(A)%A

28、= 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 （2）r_A=length(ci) r_A = 3 （3）aa=A(:,1:3)*R(1:3,4)%err=norm(A(:,4)-aa)% aa = 13 8 12 1err = 0 【例4.2-4】A=reshape(1:15,5,3);%X=null(A)%S=A*X%n=size(A,2);%l=size(X,2);%n-l=rank(A)% X = -0.4082 0.8165 -0.4082S = 1.0e-0

29、14 * 0.2665 0.2665 0.3553 0.4441 0.6217ans = 1 【例4.2-5】 （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

30、/V%A1_1=real(A1)%A2=VR*DR/VRerr1=norm(A-A1,'fro')err2=norm(A-A2,'fro') A1 = 1.0000 -3.0000 - 0.0000i 2.0000 0.6667 A1_1 = 1.0000 -3.0000 2.0000 0.6667A2 = 1.0000 -3.0000 2.0000 0.6667err1 = 4.6613e-016err2 = 4.4409e-016 4.2.3 線性方程的解 1 線性方程解的一般結(jié)論 2 除法運(yùn)算解方程【例4.2-6】（1）A=reshape(1:12,4,3

31、);%b=(13:16)'% （2）ra=rank(A)%Arab=rank(A,b)% ra = 2rab = 2 （3）xs=Ab;%xg=null(A);%c=rand(1);%ba=A*(xs+c*xg)%norm(ba-b)% Warning: Rank deficient, rank = 2, tol = 1.8757e-014.ba = 13.0000 14.0000 15.0000 16.0000ans = 1.1374e-014 3 矩陣逆【例4.2-7】（1）rng defaultA=gallery('randsvd',300,2e13,2);%x=

32、ones(300,1);%b=A*x;%cond(A)% ans = 2.0022e+013 （2）tic%xi=inv(A)*b;% ti=toc%eri=norm(x-xi)%rei=norm(A*xi-b)/norm(b)% ti = 0.2034eri = 0.0812rei = 0.0047 （3）tic;xd=Ab;%td=tocerd=norm(x-xd)red=norm(A*xd-b)/norm(b) td = 0.0125erd = 0.0274red = 7.5585e-015 4.2.4 一般代數(shù)方程的解【例 4.2-8】（1）syms tft=sin(t)2*exp(-

33、0.1*t)-0.5*abs(t);S=solve(ft,t)%<3>ftS=subs(ft,t,S)% S =0ftS =0 （2）（A）y_C=inline('sin(t).2.*exp(-0.1*t)-0.5*abs(t)','t');% （B）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)')hold off 圖 4.2-1

34、函數(shù)零點(diǎn)分布觀察圖（C）zoom on%tt,yy=ginput(5);zoom off% 圖 4.2-2 局部放大和利用鼠標(biāo)取值圖tt% tt = -2.0039 -0.5184 -0.0042 0.6052 1.6717 （D）t4,y4=fzero(y_C,0.1)%<17> t4 = 0.5993y4 = 1.1102e-016 4.3 概率分布和統(tǒng)計(jì)分析 1 二項(xiàng)分布(Binomial distribution)【例 4.3-1】N=100;p=0.5;%k=0:N;%pdf=binopdf(k,N,p);%cdf=binocdf(k,N,p);%h=plotyy(k,p

35、df,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',

36、'+','MarkerSize',4)%set(get(h(2),'Ylabel'),'String','cdf')%xlabel('k')%grid on 圖 4.3-1 二項(xiàng)分布B(100, 0.5)的概率和累計(jì)概率曲線 2 正態(tài)分布（Normal distribution）【例4.3-2】mu=3;sigma=0.5;%x=mu+sigma*-3:-1,1:3;yf=normcdf(x,mu,sigma);P=yf(4)-yf(3),yf(5)-yf(2),yf(6)-yf(1);%xd=1:

37、0.1:5;yd=normpdf(xd,mu,sigma);%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<2 text(3.8,0.6,'mu-sigma,mu+sigma') elsekk=int2str(k);text(3.8,0.6,'mu-',kk,'sigm

38、a,mu+',kk,'sigma')endtext(2.8,0.3,num2str(P(k);shgend xlabel('x');shg 圖 4.3-2 均值兩側(cè)一、二、三倍標(biāo)準(zhǔn)差之間的概率 3 各種概率分布的交互式觀察界面【例4.3-3】圖4.3-3 概率分布交互界面4.3.2 全局隨機(jī)流、隨機(jī)數(shù)組和統(tǒng)計(jì)分析 1 全局隨機(jī)流的操控表4.3-1 產(chǎn)生全局隨機(jī)流的發(fā)生器generator算 例可取字符串含義'twister'默認(rèn)的隨機(jī)數(shù)發(fā)生器例'v5uniform'MATLAB 5.0版均布隨機(jī)數(shù)發(fā)生器'v5nor

39、mal'MATLAB 5.0版正態(tài)隨機(jī)數(shù)發(fā)生例4.3-5, 例4.3-6'v4'MATLAB 4.0版隨機(jī)發(fā)生器 2 三個(gè)基本隨機(jī)數(shù)組創(chuàng)建指令【例 4.3-4】（1）rng default%<1>GRS=rand(1,25)%<2>GRS = Columns 1 through 5 0.8147 0.9058 0.1270 0.9134 0.6324 Columns 6 through 10 0.0975 0.2785 0.5469 0.9575 0.9649 Columns 11 through 15 0.1576 0.9706 0.9572

40、0.4854 0.8003 Columns 16 through 20 0.1419 0.4218 0.9157 0.7922 0.9595 Columns 21 through 25 0.6557 0.0357 0.8491 0.9340 0.6787 （2）rng default%<3>r1=randn(1,5)%<4>r2=rand(1,5)%<5>r3=randi(-3,2,1,5)%<6>r4=exprnd(0.4, 1,5)%<7>st=rng;%<8>r5=rand(1,5)%<9> r1 = 0

41、.5377 1.8339 -2.2588 0.8622 0.3188r2 = 0.0975 0.2785 0.5469 0.9575 0.9649r3 = -3 2 2 -1 1r4 = 0.7811 0.3453 0.0352 0.0932 0.0165r5 = 0.6557 0.0357 0.8491 0.9340 0.6787 （3）rng(st)%<10>rr5=rand(1,5)%<11>rr5 = 0.6557 0.0357 0.8491 0.9340 0.6787 【例 4.3-5】（1）rng(2)%N=10000;a=randn(N+2,1);A=a(

42、1:N),a(2:N+1),a(3:N+2);%A(1:4,:)%CA=corrcoef(A)% ans = -0.1242 -2.5415 0.2772 -2.5415 0.2772 -0.1960 0.2772 -0.1960 -0.1962 -0.1960 -0.1962 -0.3057CA = 1.0000 -0.0093 -0.0024 -0.0093 1.0000 -0.0093 -0.0024 -0.0093 1.0000 （2）clearrng(5)%N=10000;%A=rand(N,3);%B=randn(N,3);%C=rand(N,3);%RAB=corrcoef(A(

43、:),B(:)%RAC=corrcoef(A,C)% RAB = 1.0000 0.0017 0.0017 1.0000RAC = 1.0000 0.0083 0.0083 1.0000 （3）clearN=10000;rng(17)a=randn(1,5)%A=rand(N,3);rng(18)%b=randn(1,5)%B=rand(N,3);CAB=corrcoef(A,B)% a = -0.3951 0.1406 -1.5172 -1.8820 0.7965b = 0.2068 0.0155 0.8243 -1.6221 0.7124CAB = 1.0000 -0.0109 0.000

44、3 0.0123 -0.0045 0.0060 -0.0109 1.0000 -0.0023 0.0004 -0.0042 -0.0092 0.0003 -0.0023 1.0000 0.0247 0.0158 -0.0038 0.0123 0.0004 0.0247 1.0000 0.0167 -0.0121 -0.0045 -0.0042 0.0158 0.0167 1.0000 -0.0013 0.0060 -0.0092 -0.0038 -0.0121 -0.0013 1.0000 （4）clearN=1e4;rng default%a=rand(N,1);rng(0,'v4&

45、#39;)%b=rand(N,1);Cab=corrcoef(a,b) Cab = 1.0000 0.0031 0.0031 1.0000 表4.3-2 創(chuàng)建獨(dú)立同分布隨機(jī)序列或隨機(jī)流的不同方法匯總序號(hào)創(chuàng)建目標(biāo)創(chuàng)建特征應(yīng)用建議1接續(xù)隨機(jī)序列數(shù)組同一發(fā)生器；同一初始種子；時(shí)間上接續(xù)的內(nèi)部狀態(tài)向量最簡(jiǎn)便地產(chǎn)生滿足統(tǒng)計(jì)獨(dú)立同分布的隨機(jī)數(shù)組；用于數(shù)字信號(hào)處理中隨機(jī)的時(shí)間序列數(shù)組的產(chǎn)生2隨機(jī)數(shù)組同一發(fā)生器；同一初始種子；時(shí)間上分隔的內(nèi)部狀態(tài)向量最簡(jiǎn)便地產(chǎn)生滿足統(tǒng)計(jì)獨(dú)立同分布的隨機(jī)數(shù)組；用于Monte Carlo仿真中隨機(jī)的非時(shí)間序列數(shù)組的產(chǎn)生3多個(gè)隨機(jī)流同一發(fā)生器；不同初始種子統(tǒng)計(jì)意義上獨(dú)立性更好的隨機(jī)

46、數(shù)組；用于Monte Carlo仿真4多個(gè)隨機(jī)流不同發(fā)生器統(tǒng)計(jì)意義上獨(dú)立性最好的隨機(jī)數(shù)組；運(yùn)用于獨(dú)立同分布要求相同，而隨機(jī)流生成機(jī)理不同的訓(xùn)練集、測(cè)試集隨機(jī)數(shù)組的產(chǎn)生和Monte Carlo仿真 3 統(tǒng)計(jì)分析指令【例4.3-6】rng(0,'v5normal')%A=randn(1000,4);%AMAX=max(A)%AMIN=min(A)%CM=mean(A)%MA=mean(mean(A)%S=std(A)%var(A)-S.2%C=cov(A)diag(C)'-var(A)%p=corrcoef(A)% AMAX = 2.7316 3.2025 3.4128 3

47、.0868AMIN = -2.6442 -3.0737 -3.5027 -3.0461CM = -0.0431 0.0455 0.0177 0.0263MA = 0.0116S = 0.9435 1.0313 1.0248 0.9913ans = 1.0e-015 * 0 -0.2220 0 0C = 0.8902 -0.0528 0.0462 0.0078 -0.0528 1.0635 0.0025 0.0408 0.0462 0.0025 1.0502 -0.0150 0.0078 0.0408 -0.0150 0.9826ans = 1.0e-014 * -0.0111 -0.1554

48、-0.0888 0p = 1.0000 -0.0543 0.0478 0.0083 -0.0543 1.0000 0.0024 0.0399 0.0478 0.0024 1.0000 -0.0147 0.0083 0.0399 -0.0147 1.0000 【例4.3-7】mu=2;s=0.5;rng(22,'v5normal')%x=randn(1000,1);%<3>y=s*x+mu;%<4>z=s*(x+mu);%<5>subplot(3,1,1),histfit(x),axis(-5,5,0,100),ylabel('x')subplot(3,1,2),histfit(y),axis(-5,5,0,100),ylabel('y')subplot(3,