數(shù)學(xué)建模實驗三Lorenz模型與食餌模型

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

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

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

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

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

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

7、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求數(shù)值解的命令為：sn=NDSolve(x2-1)yx+2x*yx-Cosx=0,y0=1, y,x,0,1 畫數(shù)值解圖像的命令為： y=y/.sn1g2=Plotyx,x,0,1 比較解析解圖像與數(shù)值解圖像的命令為： Showg1,g24、求微分方程組 在初始條件,下的解，并畫出解函數(shù)的圖形。求解微分方程組的命令為：Clearx,y,txy=DSolvext+5*xt+yt

8、=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īng)出y隨x的變化關(guān)系。三、實驗準(zhǔn)備 認(rèn)真閱讀實驗?zāi)康呐c實驗材料后要正確地解讀實驗，在此基礎(chǔ)上制定實驗計劃（修改、補充或編寫程序，提出實驗思路，明確實驗步驟），為上機實驗做好準(zhǔn)備。四、實驗思路提示4.1實驗步驟1、求解問題（2）中的食餌模型的微分方程組，并畫出解的圖形和相位圖。（1）以x=800，y=100為初始值，計算x（t）,y（t），當(dāng)t0，14時的數(shù)據(jù)。繪出解的圖形，并分

9、析捕食者和被捕食者的數(shù)量變化規(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ā)現(xiàn)不能求出解析解。修改代碼如下，可以求數(shù)值解：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,0,1,RGBCo

10、lor1,0,0 捕食者和被捕食者的數(shù)量變化（2）以x為橫坐標(biāo)，y為縱坐標(biāo)繪制相位圖。根據(jù)圖形分析被捕食者數(shù)量增加(減少)對捕食者數(shù)量的影響。繪制相位圖的命令：ParametricPlotxt,yt,t,0,14 相位圖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); -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、求解以下問題（廣告的效用）： 某公司生產(chǎn)一種耐用消費品，產(chǎn)品一上市，該公司即開始做廣告，一段時