1、可達0ksksk基于MATLAB的各類混沌系統的計算機模擬混沌實驗教學平臺的設計與實現初期報告物電05級1A班 張丹偉 20050003101 摘要：本文利用數學軟件MATLAB對Lorenz系統等六個重要的混沌模型進行數值計算，同時模擬出各類混沌系統的獨特性質，如混沌吸引子，倍周期，初值敏感性，相圖，分岔圖等。通過觀察和分析上述特性，加深了我們對混沌現象的理解。關鍵詞：混沌； 微分方程； MATLAB；引言 混沌探秘混沌是非線性系統所獨有且廣泛存在的一種非周期運動形式, 其覆蓋面涉及到自然科學和社會科學的幾乎每一個分支。1972年12月29日，美國麻省理工學院教授、混沌學開創人之一E.N.洛

2、倫茲在美國科學發展學會第139次會議上發表了題為蝴蝶效應的論文，提出一個貌似荒謬的論斷：在巴西一只蝴蝶翅膀的拍打能在美國得克薩斯州產生一個龍卷風，并由此提出了天氣的不可準確預報性。為什么會出現這種情況呢？這是混沌在作怪！“混沌”譯自英語中“chaos”一詞，原意是混亂、無序，在現代非線性理論中，混沌則是泛指在確定體系中出現的貌似無規則的、類隨機的運動?；煦绗F象是普遍的，就在我們身邊，是與我們關系最密切的現象，我們就生活在混沌的海洋中。一支燃著的香煙，在平穩的氣流中緩緩升起一縷青煙，突然卷成一團團劇烈攪動的煙霧，向四方飄散；打開水龍頭，先是平穩的層流，然后水花四濺，流動變的不規則，這就是湍流；一

3、個風和日麗的夏天，突然風起云涌，來了一場暴風雨。一面旗幟在風中飄揚，一片秋葉從樹上落下，它們都在做混沌運動?？梢娀煦缡冀K圍繞在我們的周圍，一直與人類為伴?；煦绲幕靖拍?. 混沌: 目前尚無通用的嚴格的定義, 一般認為,將不是由隨機性外因引起的, 而是由確定性方程(內因)直接得到的具有隨機性的運動狀態稱為混沌。2. 相空間: 在連續動力系統中, 用一組一階微分方程描述運動, 以狀態變量(或狀態向量)為坐標軸的空間構成系統的相空間。系統的一個狀態用相空間的一個點表示, 通過該點有唯一的一條積分曲線。3. 混沌運動: 是確定性系統中局限于有限相空間的高度不穩定的運動。所謂軌道高度不穩定, 是指近鄰

4、的軌道隨時間的發展會指數地分離。由于這種不穩定性, 系統的長時間行為會顯示出某種混亂性。4. 分形和分維: 分形是 n 維空間一個點集的一種幾何性質, 該點集具有無限精細的結構, 在任何尺度下都有自相似部分和整體相似性質, 具有小于所在空間維數 n 的非整數維數。分維就是用非整數維分數維來定量地描述分形的基本性質。5. 不動點: 又稱平衡點、定態。不動點是系統狀態變量所取的一組值, 對于這些值系統不隨時間變化。在連續動力學系統中, 相空間中有一個點, 若滿足當 時, 軌跡, 則稱為不動點。6. 吸引子: 指相空間的這樣的一個點集 s (或一個子空間) , 對s鄰域的幾乎任意一點, 當時所有軌跡

5、線均趨于s, 吸引子是穩定的不動點。7. 奇異吸引子: 又稱混沌吸引子, 指相空間中具有分數維的吸引子的集合。該吸引集由永不重復自身的一系列點組成, 并且無論如何也不表現出任何周期性。混沌軌道就運行在其吸引子集中。8. 分叉和分叉點: 又稱分岔或分支。指在某個或者某組參數發生變化時, 長時間動力學運動的類型也發生變化。這個參數值(或這組參數值)稱為分叉點, 在分叉點處參數的微小變化會產生不同性質的動力學特性, 故系統在分叉點處是結構不穩定的。9. 周期解: 對于系統 , 當時,若存在 , 則稱該系統有周期解 。不動點可以看作是周期為1的解, 因為它滿足。10. 初值敏感性：對初始條件的敏感依賴

6、是混沌的基本特征，也有人用它來定義混沌：混沌系統是其終極狀態極端敏感地依賴于系統的初始狀態的系統。敏感依賴性的一個嚴重后果就在于，使得系統的長期行為變得不可預見。MATLAB中的龍格庫塔（Runge-Kutta）實現MATLAB（Matrix Laboratory）是MathWorks公司開發的，目前國際上最流行應用最廣的科學與工程計算機軟件之一。MATLAB 軟件以矩陣運算為基礎，把計算，可視化，程序設計等有機的融合在一起，具有出色的數值計算能力和強大的圖形處理功能。 基于RungeKutta法，MATLAB提供了求解微分方程數值解的函數，一般調用格式是： 其中fname 是定義的函數文件名

7、，該函數文件必須返回一個列向量。Tspan形式是t0,tf，表示求解區間，y0是初始狀態向量。這兩個函數分別采用“二階，三階RungeKutta法”和“四階，五階RungeKutta法”，并采用自適應的求解方法，即當解的變化較慢時采用較大的步長，從而使計算速度很快，當解的變化較快時步長會自動變小長，從而使計算精度很高。在MATLAB中，一般選取四階的龍格庫塔方法。Lorenz 混沌系統 美國氣象學家洛倫茲（E.N.Lorenz）于1963年在大氣科學雜志上提出第一個表現奇異吸引子的動力學系統。該混沌系統模型可以用下列微分方程組描述： 利用MATLAB數學軟件對上面微分方程求解，進行數值模擬。首

8、先建立M文件 Lorenz.m定義腳本函數，然后編程調用，其中x（1）表示x，x（2）表示y，x（3）表示z ，程序如下：function r=lorenz(t,x)global a;global b;global c;r=-c*(x(1)-x(2);a*x(1)-x(2)-x(1)*x(3);b*(x(1)*x(2)-x(3);clear;global a;global b;global c;b=8/3;c=10;t0=0,100;f0=1,1,1；for a=10:30 t,x=ode45(lorenz,t0,f0); a subplot(3,1,1); plot(t,x(:,1),r,t

9、,x(:,2),g,t,x(:,3),b); title(Lorenz 模型變量時域響應);legend(x,y,z); xlabel(t); subplot(3,1,2); plot3(x(:,1),x(:,2),x(:,3); title(Lorenz模型相圖);xlabel(x);ylabel(y);zlabel(z); grid on; subplot(3,1,3); plot(x(:,1),x(:,3); title(Lorenz模型XZ平面相圖); xlabel(x);ylabel(z); grid on; pause;end1. 固定參數b和c，設置初始值f0 和計算時間t0，通

10、過改變參數a 可以發現系統逐步進入混沌狀態的過程。2. Lorenz 吸引子 當a28時，系統已經完全進入混沌狀態，此時出現雙渦旋吸引子，如下所示：3. 倍周期：通過系數的調試可以得到Lorenz混沌的一個單倍周期和兩個多倍周期，如下：4.初值敏感性： 保持初值x0和y0不變，即x0y01，改變z0為1.001，千分之一的變化會引起系統行為的顯著改變，如下圖所示：Rossler 混沌系統Rossler 系統是化學反應系統的簡化模型，是非線性動力學中非常著名的方程，該混沌系統模型可以用下列微分方程組描述 同樣地，利用MATLAB編程求解（程序見附錄），可以對該模型進行分析。1. 逐步改變參數，觀

11、察其進入混沌狀態。Rossler吸引子： 倍周期：通過調整參數和初始值，可以得到單倍周期和2倍周期，如下圖：初值敏感性：蔡氏電路混沌系統1983年美國電學專家蔡少棠 (L. O. Chua)首次提出了著名的蔡氏電路 (Chua s circuit) ,它是混沌發展史上的一個重要里程碑. 它使人們從被動的研究混沌現象向主動的設計和控制混沌邁出了關鍵的一步.蔡氏電路是目前眾多混沌電路中最具典型代表性的一種。對于單渦旋的變形蔡氏電路的微分方程組為：當5.5，7.4， 0.25， 0.1時出現混沌狀態，如圖所示： 其初值敏感性如下所示：對于多渦旋的情況，例如可以用多項式 產生三渦卷蔡氏混沌吸引子，其無

12、量綱歸一化狀態方程可以寫為：其中，參數12.8，19.1，實驗中固定b1.1，c0.45。1. 改變參數a，觀察該電路模型的進入混沌狀態的過程當a0.6時，呈現三渦旋吸引子如下：其初值敏感性：Duffing 方程杜芬（Duffing）方程指的是非線性振子的間諧受迫振動方程： 該寫為自洽的微分方程組為： 選取參數 a2.09， b0.1，c0.5，可以得到該混沌模型的相圖：Logistic方程Logistic方程是離散映射方程：，n=0,1,2,; 其相圖如下：其分岔圖如下： Henon 方程 其中，n=0, 1, 2, ; a=1.28, b=0.3.其相圖如下： 總結與感想在這次科技制作中，

13、我主要的工作是對在已有的混沌系統方程基礎上，利用MATLAB進行編程，模擬混沌系統。由于混沌對初值和參數的依賴性很大，所以我盡可能地通過網上或者書籍尋找這些參數，減少調試時間；當然，這樣不能很好地探究整個混沌系統的特性，只能停留在理論的層面上。因此，在接下來的工作里，一方面繼續學習混沌理論（如吸引子，倍周期，Lyapunov exponent等），加深對混沌的認識；另一方面，著手混沌應用，如混沌同步控制和圖像加密通信。參考書： 1. The Essence of Chaos（混沌的本質） E.N. 洛侖茲， 氣象出版社，北京19972. 精通MATLAB 7 王正林 劉 明， 電子工業出版社，

14、北京2007 附錄資料：不需要的可以自行刪除各類濾波器的MATLAB程序理想低通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,F

15、acecolor,interp,Edgecolor,none,Facelighting,phong); 二、理想高通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(ui

16、nt8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth低通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(

17、Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth高通濾波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r(

18、i,j); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯低通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=f

19、reqspace(size(IA),meshgrid);D=100/size(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);fig

20、uresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯高通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1-exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y

21、.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 梯形低通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt

22、(f1.2+f2.2);Hd=zeros(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D1-r(i,j)/(D1-D0); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong);

23、 梯形高通濾波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=ones(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D0-r(i,j)/(D0-D1); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot

24、(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 用其他方法編寫的理想低通、理想高通、Butterworth低通、同態濾波程序理想低通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);d0=50;u0=floor(N1/2)+1;v0=floo

25、r(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d=d0 h=1; else h=0; end y(i,j)=g(i,j)*h; endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); 理想高通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);n=2;d0=10;u0=floor(N1/2)+1;v0=floor(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d=d0 h=0; else h=1; end y(