數學建模實驗三Lorenz模型與食餌模型

1、可編輯數學建模實驗三 Lorenz模型與食餌模型一、實驗目的1、學習用Mathematica求常微分方程的解析解和數值解，并進行定性分析；2、學習用MATLAB求常微分方程的解析解和數值解，并進行定性分析。二、實驗材料2.1問題圖3.3.1是著名的洛侖茲(E.N.Lorenz)混沌吸引子，洛侖茲吸引子已成為混沌理論的徽標，好比行星軌道圖代表著哥白尼、開普勒理論一樣。洛侖茲是學數學出身的，1948年起在美國麻省理工學院（MIT）作動力氣象學博士后工作，1963年他在大氣科學雜志上發(fā)表的論文確定性非周期流是混沌研究史上光輝的著作。以前科學家們不自覺地認為微分方程的解只有那么幾類：1)發(fā)散軌道；2)

2、不動點；3)極限環(huán) ；4)極限環(huán)面。除此以外，大概沒有新的運動類型了，這是人們的一種主觀猜測，誰也沒有給出證明。事實上這種想法是非常錯誤的。1963年美國麻省理工學院氣象科學家洛侖茲給出一個具體模型，就是著名的Lorenz模型，清楚地展示了一種新型運動體制：混沌運動，軌道既不收斂到極限環(huán)上也不跑掉。而今Lorenz 模型在科學與工程計算中經常運用的問題。例如，數據加密中。我們能否繪制出洛侖茲吸引子呢？圖3.3.1 洛侖茲(E.N.Lorenz)混沌吸引子假設狐貍和兔子共同生活在同一個有限區(qū)域內，有足夠多的食物供兔子享用，而狐貍僅以兔子為食物.x為兔子數量,y表狐貍數量。假定在沒有狐貍的情況下，

3、兔子增長率為400。如果沒有兔子，狐貍將被餓死，死亡率為90。狐貍與兔子相互作用的關系是，狐貍的存在使兔子受到威脅，且狐貍越多兔子增長受到阻礙越大，設增長的減小與狐貍總數成正比，比例系數為0.02。而兔子的存在又為狐貍提供食物，設狐貍在單位時間的死亡率的減少與兔子的數量成正比，設比例系數為0.001。建立數學模型，并說明這個簡單的生態(tài)系統(tǒng)是如何變化的。2.2預備知識1、求解常微分方程的Euler折線法求初值問題 （12.1）在區(qū)間上的數值解，并在區(qū)間插入了結點。由導數的定義，即微商。（右端稱為差商）從而可在每個結點上用差商來近似替代導數，將微分方程轉化為代數方程組（此處的代數方程組常稱為差分方

4、程） ,加上初值條件則可確定一組解。求解這一差分方程即可得到微分方程初值問題的數值解。變形上述方程有 ,記,，從而，則有 這就是求解微分方程初值問題的歐拉(Euler)折線法。之所以稱為歐拉折線法是因為：就幾何角度而言，所求得的近似解是初值問題精確解的折線逼近，而且此折線的起點是初值條件所對應的點。2、微分方程的Mathematica求解（1）求解命令 有兩個命令：DSolve 與NDSolve。命令格式分別為 DSolve方程，y，x NDSolve方程，y，x，xl，x2。 其中方程必須為微分方程及相應初始條件，x，xl，x2說明要給出數值解的范圍為區(qū)間x1，x2。（2）使用的注意事項 方

5、程中的函數應寫成完整形式y(tǒng)x,以表明y是x的函數； 方程應寫成=的形式； 重復使用時，應隨時清除要涉及變量的以前定義，方法是Cleary； 使用NDSolve時，所加初始條件的個數應等于微分方程的階數，同時方程中也不含其它參數，否則給不出正確結果。（3）解的表示形式Mathematica給出的微分方程的解是以純函數（或數學中的算子）定義的形式給出的，例如：DSolveyx+ 3*yx=2x,y,x的結果是3、微分方程的MATLAB求解（1）求解析解命令dsolve；（2）求數值解命令ODE或 Simulink。2.3建立模型問題（1）的洛侖茲吸引子可以用下面的微分方程得到，著名的Lorenz

6、模型的狀態(tài)方程可表示為若令 且初值為，e 為一個小常數，假設。求微分方程的數值解，并繪制出時間曲線與相空間曲線。問題（2）是著名的食餌模型，數學模型為 2.4練習題1、求解微分方程的通解。求解的Mathematica命令為： DSolveyx+2*x*yx= x*E(-x2),y,x 或者 DSolveDyx,x+2*x*yx= x*E(-x2),y,x2、求微分方程在初始條件下的特解。應給出的命令為：DSolvex*yx+ yx-Ex=0,y1=2E,y,x3、求在初始條件下的特解，并畫出解的圖形。要求分別求解析解與數值解并作比較。清除要涉及變量的命令為：Clearx,y求解析解的命令為：s

7、c=DSolve(x2-1)yx+2x*yx-Cosx=0,y0=1,y,x 畫解析解圖像的命令為：y=y/.sc1g1=Plotyx,x,0,1,PlotStyle-RGBColor1,0,0注：也可將畫圖范圍變?yōu)镻lotyx,x,0,4求數值解的命令為：sn=NDSolve(x2-1)yx+2x*yx-Cosx=0,y0=1, y,x,0,1 畫數值解圖像的命令為： y=y/.sn1g2=Plotyx,x,0,1 比較解析解圖像與數值解圖像的命令為： Showg1,g24、求微分方程組 在初始條件,下的解，并畫出解函數的圖形。求解微分方程組的命令為：Clearx,y,txy=DSolvex

8、t+5*xt+yt=Et,yt-xt-3*yt=0,x0=1,y0=0,x,y,t畫解的相位圖的命令為：y=y/.xy1;x=x/.xy1;ParametricPlotxt,yt,t,0,3,PlotRange-10,2,0,5注：圖中反應出y隨x的變化關系。三、實驗準備 認真閱讀實驗目的與實驗材料后要正確地解讀實驗，在此基礎上制定實驗計劃（修改、補充或編寫程序，提出實驗思路，明確實驗步驟），為上機實驗做好準備。四、實驗思路提示4.1實驗步驟1、求解問題（2）中的食餌模型的微分方程組，并畫出解的圖形和相位圖。（1）以x=800，y=100為初始值，計算x（t）,y（t），當t0，14時的數據。

9、繪出解的圖形，并分析捕食者和被捕食者的數量變化規(guī)律。可以先用下面的命令求解析解：Clearx,y,txy=DSolvext=4*xt-0.02*xt*yt, yt=-0.9*yt+0.001*xt*yt,x0=800, y0=100,x,y,t注：可以發(fā)現不能求出解析解。修改代碼如下，可以求數值解：Clearx,y,txy=NDSolvext=4*xt-0.02*xt*yt, yt=-0.9*yt+0.001*xt*yt,x0=800, y0=100,x,y,t,0,14繪出解的圖形：y=y/.xy1;x=x/.xy1;Plotxt,yt,t,0,14,PlotStyle-RGBColor0,

10、0,1,RGBColor1,0,0 圖3.3.2 捕食者和被捕食者的數量變化（2）以x為橫坐標，y為縱坐標繪制相位圖。根據圖形分析被捕食者數量增加(減少)對捕食者數量的影響。繪制相位圖的命令：ParametricPlotxt,yt,t,0,14 圖3.3.3 相位圖2、用MATLAB求解問題（1）中Lorenz 模型的微分方程。（1）打開MATLAB的編輯器；（2）在編輯器中用下面的幾個語句描述微分方程，并將其保存在lorenzeq.m的m文件中：f unction xdot = lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);

11、-x(1)*x(2)+28*x(2)-x(3);（3）新建命令文件：t_final=100; x0=0;0;1e-10; t,x=ode45(lorenzeq,0,t_final,x0); plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3); axis(10 40 -20 20 -20 20);繪制出時間曲線與相空間曲線，如下圖所示。圖3.3.4時間曲線與相空間曲線4.2思考問題1、運用Mathematica求解Lorenz 模型的微分方程組，從而了解系統(tǒng)狀態(tài)是如何變化的。2、求解以下問題（廣告的效用）： 某公司生產一種耐用消費品，產品一上市，該公司即開始做廣告，一段時期的