非線性動(dòng)力系統(tǒng) MatLab編程簡(jiǎn)介

1、M a t L a b 編 程 簡(jiǎn) 介一、先介紹 MatLab 里幾個(gè)挺好用,但是大家以前可能沒(méi)用到的命令,用過(guò) MatLab 的從例子都應(yīng)該能看明白 用法,所以我基本不加注釋。1、 求解初等方程:022=+solve('x2+w2'ans =i*w-i*w2、 求解微分方程:x x2-= , 以及求函數(shù)的導(dǎo)數(shù)。xt=dsolve('D2x = -w2*x','t',yt=diff(xt,'t'xt =C1*sin(w*t+C2*cos(w*tyt =C1*cos(w*t*w-C2*sin(w*t*w3、符號(hào)作圖: sin(+=t

2、 A ysyms t;%這個(gè)命令用來(lái)聲明下面要用到的符號(hào)變量。ezplot(sin(t+pi/6,0,2*pi;axis(0,2*pi,-1,1; 4、符號(hào)作帶等高線曲面圖:2222121x y z += syms x y w;ezsurfc(y*y/2+2*x*x,-4,4,-8,8; 5、畫向量場(chǎng) x y x Q y y x P sin , (, , (-=u=-4:0.3:4;x,y=meshgrid(u;%這個(gè)命令用來(lái)先生成一個(gè)數(shù)據(jù)網(wǎng)格以畫向量場(chǎng)。quiver(x,y,y,-sin(x; 6、極坐標(biāo)作圖:2-+=e r Theta=0:0.01:5*pi; polar(Theta,1.

3、/sqrt(1+exp(-2*Theta,'r' 7、畫等高線:222 1(y x x z +-=syms x y;ezcontour(x*x*(x-1*(x-1+y*y,-0.3,1.3,-0.5,0.5; MatLab 中查在線幫助。 ezmesh 畫網(wǎng)格圖ezmeshc 畫帶等高線的網(wǎng)格圖ezsurf 畫曲面圖ezsurfc 畫帶等高線的曲面圖ezplot 符號(hào)曲線圖ezplot3 三維符號(hào)曲線圖ezpolar 極坐標(biāo)曲線圖ezcontour 畫等高線圖ezcontourf 畫填充等高線圖plot 數(shù)值二維圖poar 數(shù)值極坐標(biāo)圖plot3 數(shù)值三維圖mesh 數(shù)值網(wǎng)格圖

4、meshc 數(shù)值帶等高線圖surf 數(shù)值曲面圖surfc 數(shù)值帶等高線圖quiver 矢量圖solve 符號(hào)解代數(shù)方程dsolve 符號(hào)解微分方程meshgrid 產(chǎn)生數(shù)據(jù)網(wǎng)格gradient 求數(shù)值梯度三、微分方程的數(shù)值解法微分方程 , (u t f u = ,其中 u 為矢量。以下以簡(jiǎn)諧振動(dòng)方程為例:-=x y yx 4記 (m m t u u =1、歐拉法: h u t f u u m m m m , (1+=+ (1 h=0.001;N=216;t=0:h:N*h;x=zeros(size(t;y=zeros(size(t;%可以直接取 x=zeros(1,N+1;但是一來(lái)容易范少算一

5、個(gè)點(diǎn)的錯(cuò)誤,二來(lái)下次改 t 的長(zhǎng) %度時(shí) , 還得記得改 x 的長(zhǎng)度,這樣編程 x 的長(zhǎng)度隨時(shí)間變量 t 自動(dòng)取比較好x(1=1;y(1=0;for i=1:N%這里的循環(huán)次數(shù)極容易弄錯(cuò),要么多算要么少算一個(gè)點(diǎn),這種錯(cuò)誤往往還查不出來(lái)。x(i+1=x(i+y(i*h;y(i+1=y(i-4*x(i*h;end;plot(x,y 2、龍格庫(kù)塔法: 22(6143211m m m m m m h u u +=+ (B.2 +=+=+=, ( 21, 21( 21, 21( , (3423121h u h t f h u h t f h u h t f u t f m m m m m m m m m

6、 m m m m m m (B.3 h=0.001;h2=h/2;h6=h/6;N=219;%此處設(shè)置兩個(gè)變量 h2和 h6是為了減少下面的大循環(huán)中的除法的次數(shù),記住一個(gè)原則:如果在循環(huán)里有某個(gè)量 %與循環(huán)變量無(wú)關(guān) , 那么所有關(guān)于這個(gè)量的運(yùn)算能算好的都先算好。t=0:h:N*h;x=zeros(size(t;y=zeros(size(t;x(1=1;y(1=0;for i=1:Nux=x(i; uy=y(i; wx1=uy; wy1=-4*ux;ux=x(i+wx1*h2; uy=y(i+wy1*h2; wx2=uy; wy2=-4*ux;ux=x(i+wx2*h2; uy=y(i+wy2*

7、h2; wx3=uy; wy3=-4*ux;ux=x(i+wx3*h; uy=y(i+wy3*h; wx4=uy; wy4=-4*ux;x(i+1=x(i+(wx1+2*wx2+2*wx3+wx4*h6;y(i+1=y(i+(wy1+2*wy2+2*wy3+wy4*h6;end;%在上面的循環(huán)里,我設(shè)了兩個(gè)臨時(shí)變量 ux , uy ,目前大家能看到的是兩個(gè)優(yōu)點(diǎn),一個(gè)是程序結(jié)構(gòu)非常清楚, %另一個(gè)優(yōu)點(diǎn)是這個(gè)程序的可移植性非常強(qiáng),如果算另一個(gè)方程,只要作極小的改動(dòng)就行了,見下一個(gè)例子 %至于用兩個(gè)臨時(shí)變量是否能加快運(yùn)行速度,這個(gè)例子體現(xiàn)不出來(lái),而下一個(gè)例子,則會(huì)有顯著區(qū)別plot(x,y N=21

8、6個(gè)點(diǎn),圖像已 經(jīng)成了橢圓環(huán)了,而后面的龍格-庫(kù)塔法,算了 N=219,是前者的 8倍,圖像仍然是一個(gè)漂亮的橢圓3、程序的優(yōu)化下面我們用龍格-庫(kù)塔法來(lái)畫下面方程的軌線:=+=y x y x y x cos cos sin (1 不良編程的典范tich=0.001;N=219;t=0:h:N*h;x=zeros(size(t;y=zeros(size(t;x(1=1;y(1=0;for i=1:Nwx1=sin(y(i+cos(x(i; wy1=x(i*cos(y(i;wx2=sin(y(i+wy1*h/2+cos(x(i+wx1*h/2; wy2=(x(i+wx1*h/2*cos(y(i+wy

9、1*h/2; wx3=sin(y(i+wy2*h/2+cos(x(i+wx2*h/2; wy3=(x(i+wx2*h/2*cos(y(i+wy2*h/2; wx4=sin(y(i+wy3*h+cos(x(i+wx3*h; wy4=(x(i+wx3*h*cos(y(i+wy3*h; x(i+1=x(i+(wx1+2*wx2+2*wx3+wx4*h/6;y(i+1=y(i+(wy1+2*wy2+2*wy3+wy4*h/6;end;plot(x,ytoc (2 常量的使用tich=0.001;h2=h/2;h6=h/6;N=219;t=0:h:N*h;x=zeros(size(t;y=zeros(s

10、ize(t;x(1=1;y(1=0;for i=1:Nwx1=sin(y(i+cos(x(i; wy1=x(i*cos(y(i;wx2=sin(y(i+wy1*h2+cos(x(i+wx1*h2; wy2=(x(i+wx1*h2*cos(y(i+wy1*h2; wx3=sin(y(i+wy2*h2+cos(x(i+wx2*h2; wy3=(x(i+wx2*h2*cos(y(i+wy2*h2; wx4=sin(y(i+wy3*h+cos(x(i+wx3*h; wy4=(x(i+wx3*h*cos(y(i+wy3*h; x(i+1=x(i+(wx1+2*wx2+2*wx3+wx4*h6;y(i+1

11、=y(i+(wy1+2*wy2+2*wy3+wy4*h6;end;plot(x,ytocElapsed time is 13.559000 seconds. 現(xiàn)在我們看到,常量 h2和 h6的引入使得我們的一條軌道的計(jì)算時(shí)間減少了 1.5秒。我們還可以進(jìn)一步改進(jìn),引進(jìn)局部變量。(3使用局部變量tich=0.001;h2=h/2;h6=h/6;N=219;t=0:h:N*h;x=zeros(size(t;y=zeros(size(t;x(1=1;y(1=0;for i=1:Nux=x(i; uy=y(i; wx1=sin(uy+cos(ux; wy1=ux*cos(uy; ux=x(i+wx1*

12、h2; uy=y(i+wy1*h2; wx2=sin(uy+cos(ux; wy2=ux*cos(uy; ux=x(i+wx2*h2; uy=y(i+wy2*h2; wx3=sin(uy+cos(ux; wy3=ux*cos(uy; ux=x(i+wx3*h; uy=y(i+wy3*h; wx4=sin(uy+cos(ux; wy4=ux*cos(uy; x(i+1=x(i+(wx1+2*wx2+2*wx3+wx4*h6;y(i+1=y(i+(wy1+2*wy2+2*wy3+wy4*h6;end;plot(x,ytocElapsed time is 11.958000 seconds. (i時(shí)

13、間上的節(jié)省:我們看到,時(shí)間又節(jié)省了 1.6秒,我們前后共節(jié)約時(shí)間 3.1秒,不要小看這 3.1秒,這只是 計(jì)算一條軌道,如果我們要畫出一個(gè)參數(shù)為 1, 0; 1, 0的舌頭,如果取步長(zhǎng)為 0.001(這個(gè)不算很精確,且不 管計(jì)算 Poincare 截面所增加的大量時(shí)間,就算一個(gè)網(wǎng)格點(diǎn)算一條軌道,那么就要增加 1000×1000×3.1/86400=35.8796(天! (ii可移植性:我們看到引入臨時(shí)變量后,后面的這個(gè)方程和前面的例子改動(dòng)的地方很少,而且極有規(guī)律 (iii 這個(gè)例子中,臨時(shí)變量的引入,之所以能使計(jì)算的時(shí)間縮短,在于龍格-庫(kù)塔法中,現(xiàn)在的例子里每個(gè) ux 和

14、uy 都要使用兩次,因此還有一個(gè)原則:如果在循環(huán)體中某個(gè)與循環(huán)變量有關(guān)的量會(huì)用到不止一次,那 就應(yīng)該增加一個(gè)臨時(shí)變量。道理弄明白了,那就只看大家的發(fā)揮了。四、一個(gè)具體的系統(tǒng)下面講一個(gè)具體的系統(tǒng):強(qiáng)迫 Duffing 方程,其中的程序由上面的討論可以知道,可移植性很好,大家 可以自由使用,但請(qǐng)不要改動(dòng)其中的版權(quán)說(shuō)明部分 (一般發(fā)布源程序的人,都要寫上這句話的! ,程序因?yàn)?目前是在 Word 里運(yùn)行,速度比較慢,所以運(yùn)行長(zhǎng)度都取得比較小,大家最好直接在 MatLab 里設(shè)置大一點(diǎn)的 計(jì)算長(zhǎng)度調(diào)用執(zhí)行。這里的幾個(gè)程序都是寫成函數(shù)直接調(diào)用的,在 MatLab 里面的調(diào)用格式,和在這里用的 是一樣的。

15、另外我不保證程序的正確性,那是你們的事情,切記切記!t F x x x x cos 3=+ 或 +-=t F x x y y y x cos 31、 軌線圖:Duffing(x0,y0,a,b,c,F,w,N,Dirx0,y0顯然為初始值, a,b,c,w 分別對(duì)應(yīng)于參數(shù) , , , , N 為計(jì)算點(diǎn)數(shù), Dir 為方向。第 11頁(yè)(共 16頁(yè)a=0.3;b=1;c=1;w=1.2;N=217;F=0.2; subplot(2,2,1;Duffing(-1.1,0,a,b,c,F,w,N;title('F=0.2'F=0.27;subplot(2,2,2;Duffing(-1.

16、1,0,a,b,c,F,w,N;title('F=0.27'F=0.29;subplot(2,2,3;Duffing(-1.1,0,a,b,c,F,w,N;title('F=0.29'F=0.32;subplot(2,2,4;Duffing(-1.1,0,a,b,c,F,w,N;title('F=0.32' 2、 變步長(zhǎng)龍格庫(kù)塔法軌線圖: BBCDuffing(x0,y0,a,b,c,F,w,Nx0,y0顯然為初始值, a,b,c,w 分別對(duì)應(yīng)于參數(shù) , , , , N 為次數(shù)。a=0.3;b=1;c=1;w=1.2;N=217;F=0.2; s

17、ubplot(2,2,1;BBCDuffing(-1.1,0,a,b,c,F,w,N;title('F=0.2'F=0.27;subplot(2,2,2;BBCDuffing(-1.1,0,a,b,c,F,w,N;title('F=0.27'F=0.29;subplot(2,2,3;BBCDuffing(-1.1,0,a,b,c,F,w,N;title('F=0.29'F=0.32;subplot(2,2,4;BBCDuffing(-1.1,0,a,b,c,F,w,N;title('F=0.32' 變步長(zhǎng)的效果從前三個(gè)圖中看不出來(lái)

18、,因?yàn)?前三個(gè)圖是周期的,而第四個(gè)圖就有明顯的差別了,我們作了 同樣的迭代次數(shù),但是 BBC 的軌道顯然走了更長(zhǎng)的距離。3、 Poincare 映射:DuffingFlash(x0,y0,a,b,c,F,w,Delta,Count其它參數(shù)同上, Delta 為頻閃步長(zhǎng),當(dāng) Delta 等于驅(qū)動(dòng)外力周期時(shí),就是 Poincare 截面, Count 為頻閃取點(diǎn) 數(shù),計(jì)算步長(zhǎng)是程序自己內(nèi)部自適應(yīng)調(diào)整的,無(wú)需設(shè)置。a=0.3;b=1;c=1;w=1.2;Count=100;F=0.2; Delta=2*pi/w;subplot(2,2,1;DuffingFlash(-1.1,0,a,b,c,F,w,

19、Delta,Count; F=0.27; Delta=2*pi/w;subplot(2,2,2;DuffingFlash(-1.1,0,a,b,c,F,w,Delta,Count; F=0.29; Delta=2*pi/w;subplot(2,2,3;DuffingFlash(-1.1,0,a,b,c,F,w,Delta,Count; F=0.32; Delta=2*pi/w;subplot(2,2,4;DuffingFlash(-1.1,0,a,b,c,F,w,Delta,500;第 12頁(yè)(共 16頁(yè) 為確定前面三種情況確實(shí)是周期函數(shù),上面的 Poincare 映射圖給了我們清楚的圖像。注

20、意第一個(gè)圖其實(shí)只有 一個(gè)點(diǎn)。3、 功率譜:DuffingGLP(x0,y0,a,b,c,F,w,N,CountN 是計(jì)算功率譜的長(zhǎng)度,盡可能設(shè)置成 2的冪,這是由快速傅立葉變換決定的。 Count 為要顯 示的功率譜的點(diǎn)數(shù)。程序是調(diào)用快速傅立葉變換,其實(shí) MatLab 里有專門的計(jì)算功率譜的函數(shù) PSD ,還有些別的函數(shù),因?yàn)樽罱鼪](méi)有時(shí)間去試,有時(shí)間的人自己去看 MatLab 的在線幫助。a=0.3;b=1;c=1;w=1.2;N=218;Count=100;F=0.2; subplot(2,2,1;DuffingGLP(-1.1,0,a,b,c,F,w,N,Count;F=0.27;subp

21、lot(2,2,2;DuffingGLP(-1.1,0,a,b,c,F,w,N,Count;F=0.29;subplot(2,2,3;DuffingGLP(-1.1,0,a,b,c,F,w,N,Count;F=0.32;subplot(2,2,4;DuffingGLP(-1.1,0,a,b,c,F,w,N,Count;第 13頁(yè)(共 16頁(yè) 功率譜圖表明,解不斷出現(xiàn)新的分頻部分,最后呈現(xiàn)沒(méi)有周期的運(yùn)動(dòng)。系統(tǒng)從什么時(shí)候開始軌線變得不是周期的呢?這就用到我們的分叉圖。4、 分叉圖:DuffingFCH(x0,y0,a,b,c,w,FMin,FMax,dF,CountFMin , FMax , dF

22、 分別為可調(diào)整參數(shù)的最小值,最大值和計(jì)算步長(zhǎng)。這三個(gè)值決定了橫向所要 顯示的點(diǎn)數(shù), Count 為縱向所要顯示的點(diǎn)數(shù),也就是 Poincare 截面的點(diǎn)數(shù)。目前下面的圖是有問(wèn) 題的, F<0.2時(shí),有好幾條線,那是因?yàn)榍懊鎸懗绦虻臅r(shí)候,暫態(tài)過(guò)程忽略過(guò)少,現(xiàn)在程序內(nèi)部 已經(jīng)增加了暫態(tài)過(guò)程忽略的長(zhǎng)度,但是這個(gè)程序運(yùn)行太費(fèi)時(shí)間,所以請(qǐng)大家自己去運(yùn)行,我就 不再運(yùn)行一次了。a=0.3;b=1;c=1;w=1.2;FMin=0;FMax=0.4;dF=0.001;Count=100;DuffingFCH(-1.2,0,a,b,c,w,FMin,FMax,dF,Count;第 14頁(yè)(共 16頁(yè)第 15頁(yè)(共 16頁(yè) 從圖中可以看到,系統(tǒng)大概在 3. 0 F 的地方開始變得沒(méi)有周期了。5、 Liapunov 指數(shù) DuffingLZS(x0,y0,a,b,c,w,F,Count;因?yàn)槌绦虍?dāng)時(shí)寫的時(shí)候只為教學(xué)用的,所以除了系統(tǒng)的幾個(gè)參數(shù)