數(shù)值分析 Lorenz問題與混沌

1、991790 力92 王為實(shí)驗(yàn)6.1（Lorenz問題與混沌）問題提出：考慮著名的Lorenz方程其中s,r,b為變化區(qū)域有一定限制的實(shí)參數(shù)，該方程形式簡(jiǎn)單，表面上看并無(wú)驚人之處，但由該方程揭示出的許多現(xiàn)象，促使“混沌”成為數(shù)學(xué)研究的嶄新領(lǐng)域，在實(shí)際應(yīng)用中也產(chǎn)生了巨大的影響。實(shí)驗(yàn)要求：（1）請(qǐng)讀者找出Lorenz方程與上述程序中使用的方程間的關(guān)系。（2）對(duì)目前取定的參數(shù)值SIGMA、RHO和BETA，選取不同的初始值y0（當(dāng)前程序中的y0是在坐標(biāo)原點(diǎn)），運(yùn)行上述的程序，觀察計(jì)算的結(jié)果有什么特點(diǎn)？解的曲線是否有界？解的曲線是不是周期的或趨于某個(gè)固定的點(diǎn)？（3）在問題允許的范圍內(nèi)適當(dāng)改變其中的參數(shù)

2、值SIGMA、RHO和BETA，再選取不同的初始值y0，運(yùn)行上述的程序，觀察并記錄計(jì)算的結(jié)果有什么特點(diǎn)？是否發(fā)現(xiàn)什么不同的現(xiàn)象。程序清單：主程序：%BLZ Plot the orbit around the Lorenz chaotic attractorclfclc%Solve the ordinary differential equation describing the %Lorenz chaotic attractor.The equation are difined in%an M-file,blzeq.m%The value of the global parameters ar

3、eglobal SIGMA RHO BETASIGMA=10.;RHO=10.;BETA=20.;%The graphics axis limits are set to values known to%contain the solution.axis(0 50 -30 30 -30 30)view(3)hold ontitle('Lorenz Attractor')y0=0 0 eps;tfinal=100;t,y=ode23('blzeq',0 tfinal,y0);plot3(y(:,1),y(:,2),y(:,3)blzeq函數(shù)程序%BLZEQ Equ

4、ation of the Lorenz chaotic attractor%ydot=lorenzeq(t,y)%The differential equation is written in almost linear form.function ydot=blzeq(t,y);global SIGMA RHO BETAA=-BETA 0 y(2) 0 -SIGMA SIGMA -y(2) RHO -1;ydot=A*y;實(shí)驗(yàn)結(jié)果及其分析：題中的方程與程序中的方程的關(guān)系是變量進(jìn)行了輪換，x換成了y，y換成了z，z換成了x。原點(diǎn)為原方程的一個(gè)奇點(diǎn)，當(dāng)初始位置稍稍偏離原點(diǎn)如取為0，eps，0，（

5、按原方程中的順序，下同）得到的圖像如下：分析：這是一個(gè)典型的奇怪吸引子的圖像，曲線有界，但他不收斂于某一點(diǎn)也不是周期的，而是在兩個(gè)位置附近來回的跳躍。取初始位置分別為10，10，10，50，50，50得到的圖像如下：分析：個(gè)初始變量值相同時(shí)，曲線總是被吸引回奇怪吸引子附近作來回跳躍，在初始變量值取道200，200，200，-200，-200，-200時(shí)，依然如此。下面分別考察初始值的每個(gè)分量變化對(duì)圖像的影響：y分量：0，3，0 0，7，00，7.69516，0 0，7.69517，00，14.23，0 0，14.24，0分析：從上面可以看見，隨著初始x值的增大，奇怪吸引子中曲線在其附近來回跳躍

6、的兩個(gè)位置中的一個(gè)吸引力變?nèi)?，另一個(gè)吸引力變強(qiáng)，然后在初始x取某一特定值時(shí)，一個(gè)位置喪失吸引力，另一位置則將曲線完全吸引過來變成普通吸引子。初始x繼續(xù)增大到某一特定值，情況又會(huì)變回來。對(duì)初始x值單獨(dú)變化的情況也有類似的現(xiàn)象。這所明在空間存在一些區(qū)域，當(dāng)初始位置位于這些區(qū)域外時(shí)解將出現(xiàn)奇怪吸引子的性質(zhì)，而在這些區(qū)域以內(nèi)解將呈現(xiàn)普通吸引子的性質(zhì)。Z分量：0，0，20分析：從圖上可以看出解的曲線為一直線，這可以從方程的角度來解釋。當(dāng)x=0,y=0時(shí)在方程中dx/dt=0,dy/dt=0，x，y 方向的值不發(fā)生變化，僅z方向的值變化，因此解為一直線。當(dāng)SIGMA=10，RHO=20，BETA=10，初

7、始值取0，eps，0時(shí)，圖像如下：對(duì)初值進(jìn)行調(diào)整沒有發(fā)現(xiàn)奇怪吸引子的出現(xiàn)。只調(diào)整BETA變量情況如下：SIGMA=10，RHO=28，BETA=9.6/3 0，eps，0SIGMA=10，RHO=28，BETA=9.7/3 0，eps，0SIGMA=10，RHO=28，BETA=15/3 0，eps，0改變SIGMA、RHO的值也有類似的現(xiàn)象。實(shí)驗(yàn)6.2（剛性問題）問題提出：考慮下面的初值問題，其中實(shí)驗(yàn)內(nèi)容：剛性比是衡量問題困難程度的重要指標(biāo)，針對(duì)問題合理選擇求解剛性問題的方法很重要，Matlab中提供了豐富的函數(shù)求解剛性方程，請(qǐng)讀者嘗試不同方法求解上述方程。實(shí)驗(yàn)要求：（1）用Runge-Ku

8、tta算法求解方程。（2）分別用剛性問題的算法和一般問題的算法求解方程，與（1）比較他們的計(jì)算結(jié)果和計(jì)算時(shí)間，并分析它們的精度。（3）在t=0.1,0.5,1,2,4,8,10等處，計(jì)算解的剛性比。（4）嘗試編寫隱式Runge-Kutta算法和非線性方法的程序，計(jì)算方程的解并與前面的計(jì)算結(jié)果相比較。程序清單：用龍格庫(kù)塔法求解方程的程序如下：主程序：clfclcaxis(-0.5 1.5 0.8 2.2, -0.4 0.4)view(3)hold ony0=1 1 0;tfinal=10;tic;t,y=ode23('func1',0 tfinal,y0);tocplot3(y(

9、:,1),y(:,2),y(:,3)函數(shù)func1.m:function ydot=func1(t,y);global SIGMA RHO BETAA=-0.013 -1000*y(1) 0 0 -2500*y(3) 0 -0.013 -1000*y(1) -2500*y(2);ydot=A*y;用Matlab中專門處理剛性問題的算法求解方程的程序如下：clfclcaxis(-0.5 1.5 0.8 2.2, -0.4 0.4)view(3)hold ony0=1 1 eps;tfinal=10;tic;t,y=ode23s('func1',0 tfinal,y0)tocplo

10、t3(y(:,1),y(:,2),y(:,3)自編的隱式龍格-庫(kù)塔程序如下：%一階隱式龍格-庫(kù)塔法clfclcaxis(-0.5 1.5 0.8 2.2 -0.4 0.4)view(3)hold onh=0.0001;y(1,:)=1 1 0;for i=1:1/h; z1=y(i,:);%yn z2=z1;%設(shè)置z2為yn+1的迭代初值 F=z1(1)-z2(1)-h*(0.0065*z1(1)+0.0065*z2(1)+250*(z1(1)+z2(1)*(z1(2)+z2(2) z1(2)-z2(2)-625*h*(z1(2)+z2(2)*(z1(3)+z2(3) z1(3)-z2(3)-

11、h*(0.0065*(z1(1)+z2(1)+250*(z1(1)+z2(1)*(z1(2)+z2(2)+625*(z1(2)+z2(2)*(z1(3)+z2(3); F1=-1-0.0065*h-250*h*(z1(2)+z2(2) -250*h*(z1(1)+z2(1) 0 0 -1-625*h*(z1(3)+z2(3) -625*h*(z1(2)+z2(2) -0.0065*h-250*h*(z1(2)+z2(2) -250*h*(z1(1)+z2(1)-625*h*(z1(3)+z2(3) -625*h*(z1(2)+z2(2)-1; z3=z2-(inv(F1)*F)'%迭代

12、求解 while norm(z3-z2)>10; z2=z3; F=z1(1)-z2(1)-h*(0.0065*z1(1)+0.0065*z2(1)+250*(z1(1)+z2(1)*(z1(2)+z2(2) z1(2)-z2(2)-625*h*(z1(2)+z2(2)*(z1(3)+z2(3) z1(3)-z2(3)-h*(0.0065*(z1(1)+z2(1)+250*(z1(1)+z2(1)*(z1(2)+z2(2)+625*(z1(2)+z2(2)*(z1(3)+z2(3); F1=-1-0.0065*h-250*h*(z1(2)+z2(2) -250*h*(z1(1)+z2(1

13、) 0 0 -1-625*h*(z1(3)+z2(3) -625*h*(z1(2)+z2(2) -0.0065*h-250*h*(z1(2)+z2(2) -250*h*(z1(1)+z2(1)-625*h*(z1(3)+z2(3) -625*h*(z1(2)+z2(2)-1; z3=z2-(inv(F1)*F)'end; y(i+1,:)=z3;end;plot3(y(:,1),y(:,2),y(:,3);tfinal=10;y0=y(1,:)自編的龍格-庫(kù)塔型非線性法的程序如下：%龍格-庫(kù)塔型非線性法clfclcaxis(-0.5 1.5 0.8 2.2 -0.4 0.4)view(

14、3)hold onh=0.0001;y(1,:)=1 1 0;for i=1:1/h; F=-0.013*y(i,1)-1000*y(i,1)*y(i,2) -2500*y(i,2)*y(i,3) -0.013*y(i,1)-1000*y(i,1)*y(i,2)-2500*y(i,2)*y(i,3);%求導(dǎo)數(shù)值x=y(i,1)+0.5*h*y(i,1)*F(1)/(y(i,1)-0.5*h*F(1) y(i,2)+0.5*h*y(i,2)*F(1)/(y(i,1)-0.5*h*F(2) y(i,3)+0.5*h*y(i,3)*F(3)/(y(i,3)-0.5*h*F(3);F1=-0.013*

15、x(1)-1000*x(1)*x(2) -2500*x(2)*x(3) -0.013*x(1)-1000*x(1)*x(2)-2500*x(2)*x(3);y(i+1,1)=y(i,1)+h*F1(1);y(i+1,2)=y(i,2)+h*F1(2);y(i+1,3)=y(i,3)+h*F1(3);end;plot3(y(:,1),y(:,2),y(:,3);y0=y(1,:)tic;t,y=ode23s('func1',0 10,y0)tocplot3(y(:,1),y(:,2),y(:,3),'r')實(shí)驗(yàn)結(jié)果及其分析：（1） 用matlab中的ode23求解的結(jié)果如下：占用的CPU時(shí)間為：6.9680s（2）用matlab中求解剛性問題的函數(shù)ode23s求解的結(jié)果如下（將它與ode23算得的結(jié)果畫在一起了）占用的CPU時(shí)間