非線性動力系統(tǒng) MatLab編程簡介

1、M a t L a b 編 程 簡 介一、先介紹 MatLab 里幾個挺好用,但是大家以前可能沒用到的命令,用過 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、符號作圖: sin(+=t

2、 A ysyms t;%這個命令用來聲明下面要用到的符號變量。ezplot(sin(t+pi/6,0,2*pi;axis(0,2*pi,-1,1; 4、符號作帶等高線曲面圖:2222121x y z += syms x y w;ezsurfc(y*y/2+2*x*x,-4,4,-8,8; 5、畫向量場 x y x Q y y x P sin , (, , (-=u=-4:0.3:4;x,y=meshgrid(u;%這個命令用來先生成一個數(shù)據(jù)網(wǎ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 符號曲線圖ezplot3 三維符號曲線圖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 符號解代數(shù)方程dsolve 符號解微分方程meshgrid 產(chǎn)生數(shù)據(jù)網(wǎng)格gradient 求數(shù)值梯度三、微分方程的數(shù)值解法微分方程 , (u t f u = ,其中 u 為矢量。以下以簡諧振動方程為例:-=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;但是一來容易范少算一

5、個點(diǎn)的錯誤,二來下次改 t 的長 %度時 , 還得記得改 x 的長度,這樣編程 x 的長度隨時間變量 t 自動取比較好x(1=1;y(1=0;for i=1:N%這里的循環(huán)次數(shù)極容易弄錯,要么多算要么少算一個點(diǎn),這種錯誤往往還查不出來。x(i+1=x(i+y(i*h;y(i+1=y(i-4*x(i*h;end;plot(x,y 2、龍格庫塔法: 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è)置兩個變量 h2和 h6是為了減少下面的大循環(huán)中的除法的次數(shù),記住一個原則:如果在循環(huán)里有某個量 %與循環(huán)變量無關(guān) , 那么所有關(guān)于這個量的運(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è)了兩個臨時變量 ux , uy ,目前大家能看到的是兩個優(yōu)點(diǎn),一個是程序結(jié)構(gòu)非常清楚, %另一個優(yōu)點(diǎn)是這個程序的可移植性非常強(qiáng),如果算另一個方程,只要作極小的改動就行了,見下一個例子 %至于用兩個臨時變量是否能加快運(yùn)行速度,這個例子體現(xiàn)不出來,而下一個例子,則會有顯著區(qū)別plot(x,y N=21

8、6個點(diǎn),圖像已 經(jīng)成了橢圓環(huán)了,而后面的龍格-庫塔法,算了 N=219,是前者的 8倍,圖像仍然是一個漂亮的橢圓3、程序的優(yōu)化下面我們用龍格-庫塔法來畫下面方程的軌線:=+=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的引入使得我們的一條軌道的計算時間減少了 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時

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

14、uy 都要使用兩次,因此還有一個原則:如果在循環(huán)體中某個與循環(huán)變量有關(guān)的量會用到不止一次,那 就應(yīng)該增加一個臨時變量。道理弄明白了,那就只看大家的發(fā)揮了。四、一個具體的系統(tǒng)下面講一個具體的系統(tǒng):強(qiáng)迫 Duffing 方程,其中的程序由上面的討論可以知道,可移植性很好,大家 可以自由使用,但請不要改動其中的版權(quán)說明部分 (一般發(fā)布源程序的人,都要寫上這句話的! ,程序因為 目前是在 Word 里運(yùn)行,速度比較慢,所以運(yùn)行長度都取得比較小,大家最好直接在 MatLab 里設(shè)置大一點(diǎn)的 計算長度調(diào)用執(zhí)行。這里的幾個程序都是寫成函數(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 分別對應(yīng)于參數(shù) , , , , N 為計算點(diǎn)數(shù), Dir 為方向。第 11頁(共 16頁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、 變步長龍格庫塔法軌線圖: BBCDuffing(x0,y0,a,b,c,F,w,Nx0,y0顯然為初始值, a,b,c,w 分別對應(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' 變步長的效果從前三個圖中看不出來

18、,因為,前三個圖是周期的,而第四個圖就有明顯的差別了,我們作了 同樣的迭代次數(shù),但是 BBC 的軌道顯然走了更長的距離。3、 Poincare 映射:DuffingFlash(x0,y0,a,b,c,F,w,Delta,Count其它參數(shù)同上, Delta 為頻閃步長,當(dāng) Delta 等于驅(qū)動外力周期時,就是 Poincare 截面, Count 為頻閃取點(diǎn) 數(shù),計算步長是程序自己內(nèi)部自適應(yīng)調(diào)整的,無需設(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頁(共 16頁 為確定前面三種情況確實(shí)是周期函數(shù),上面的 Poincare 映射圖給了我們清楚的圖像。注

20、意第一個圖其實(shí)只有 一個點(diǎn)。3、 功率譜:DuffingGLP(x0,y0,a,b,c,F,w,N,CountN 是計算功率譜的長度,盡可能設(shè)置成 2的冪,這是由快速傅立葉變換決定的。 Count 為要顯 示的功率譜的點(diǎn)數(shù)。程序是調(diào)用快速傅立葉變換,其實(shí) MatLab 里有專門的計算功率譜的函數(shù) PSD ,還有些別的函數(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頁(共 16頁 功率譜圖表明,解不斷出現(xiàn)新的分頻部分,最后呈現(xiàn)沒有周期的運(yùn)動。系統(tǒng)從什么時候開始軌線變得不是周期的呢?這就用到我們的分叉圖。4、 分叉圖:DuffingFCH(x0,y0,a,b,c,w,FMin,FMax,dF,CountFMin , FMax , dF

22、 分別為可調(diào)整參數(shù)的最小值,最大值和計算步長。這三個值決定了橫向所要 顯示的點(diǎn)數(shù), Count 為縱向所要顯示的點(diǎn)數(shù),也就是 Poincare 截面的點(diǎn)數(shù)。目前下面的圖是有問 題的, F<0.2時,有好幾條線,那是因為前面寫程序的時候,暫態(tài)過程忽略過少,現(xiàn)在程序內(nèi)部 已經(jīng)增加了暫態(tài)過程忽略的長度,但是這個程序運(yùn)行太費(fèi)時間,所以請大家自己去運(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頁(共 16頁第 15頁(共 16頁 從圖中可以看到,系統(tǒng)大概在 3. 0 F 的地方開始變得沒有周期了。5、 Liapunov 指數(shù) DuffingLZS(x0,y0,a,b,c,w,F,Count;因為程序當(dāng)時寫的時候只為教學(xué)用的,所以除了系統(tǒng)的幾個參數(shù)