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

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

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

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

4、程y f (x ,y)轉(zhuǎn)化為代數(shù)方程組(此處的代數(shù)方程組常稱為差分方程)y(Xk h) y(xjf(xk , y(xj) ,k 0,1, n 1 h加上初值條件則可確定一組解。求解這一差分方程即可得到微分方程初值問題的數(shù)值解。變形上述方程有y(xkh)y(xjhf 的，y(xj) ,k 0,1, n 1記xk1xk h, y(xk)yk,從而 y(xkh) yk 1 ,則有y0 y(x0),xk 1 xk h , k 0,1, n 1yk 1ykhf (xk , yk),這就是求解微分方程初值問題的歐拉(Euler)折線法。之所以稱為歐拉折線法是因為：就幾何角度而言，所求得的近似解是初值問題精

5、確解的折線逼近，而且此折線的起點是初值條件所對應(yīng)的點。2、微分方程的 Mathematica求解 (1)求解命令有兩個命令：DSolve與NDSolve。命令格式分別為DSolve方程，y, xNDSolve 方程，y, x, xl, x2。其中方程必須為微分方程及相應(yīng)初始條件， x, xl, x2說明要給出數(shù)值解的范圍為區(qū)間 x1 , x2 。(2)使用的注意事項方程中的函數(shù)應(yīng)寫成完整形式y(tǒng)x,以表明y是x的函數(shù)；方程應(yīng)寫成=的形式；重復(fù)使用時，應(yīng)隨時清除要涉及變量的以前定義，方法是 Cleary;使用NDSolve時，所加初始條件的個數(shù)應(yīng)等于微分方程的階數(shù)，同時方程中也不含其它參數(shù)，否則給

6、不出正確結(jié)果。 (3)解的表示形式Mathematica 給出的微分方程的解是以純函數(shù)(或數(shù)學(xué)中的算子)定義的形式給出的，例如： DSolveyx+ 3*yx=2x,y,x的結(jié)果是Function 依. W 上十：)+ L C 1 93：3、微分方程的 MATLAB求解 (1)求解析解命令 dsolve ； (2)求數(shù)值解命令 ODE或 Simulinko 建立模型問題(1)的洛侖茲吸引子可以用下面的微分方程得到，著名的Lorenz模型的狀態(tài)方程可表示為Xi（t） Xi（t） X2（t）X3（t）X2（t）X2（t） X3（t）X3（t）Xi（t）X2（t）X2（t） X3（t）若令10,28

7、,8/3,且初值為Xi（0）X2（0）0, X3（0）,為一個小常數(shù)，假設(shè)10 10。求微分方程的數(shù)值解，并繪制出時間曲線與相空間曲線。問題（2）是著名的食餌模型，數(shù)學(xué)模型為x 4x 0.02Xyy 0.9y 0.001xy練習(xí)題21、求解微分方程y 2xy xe X的通解。求解的Mathematica命令為：DSolveyx+2*x*yx= x*EA（-xA2）,y,x或者DSolveDyx,x+2*x*yx= x*EA（-xA2）,y,x2、求微分方程xy y eX0在初始條件yx12e下的特解。應(yīng)給出的命令為：DSolvex*yx+ yx-EAx=0,y1=2E,y,x3、求（x21）d

8、y 2xy cosx0在初始條件y（0） 1下的特解，并畫出解的圖形。要求分dx別求解析解與數(shù)值解并作比較。清除要涉及變量的命令為：Clearx,y求解析解的命令為：sc=DSolve（xA2-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（xA2-1）yx+2x*yx卜Cosx=0,y0=1, y,x,0,1畫數(shù)值解圖像的命令為：y=y/.sn1g2=Plotyx,x,0,1比較解析解圖像

9、與數(shù)值解圖像的命令為：Showg1,g24、求微分方程組dxdt dy dt5x3y 0在初始條件x(0) 1,y(0)0下的解，并畫出解函數(shù)y y(x)的圖形求解微分方程組的命令為：Clearx,y,txy=DSolvext+5*xt+yt=EAt,yt-xt-3*yt=0,x0=1,y0=0,x,y,t畫解的相位圖的命令為：y=y/.xyi;x=x/.xy1;ParametricPlotxt,yt,t,0,3,PlotRange-10,2,0,5注：圖中反應(yīng)出y隨x的變化關(guān)系。三、實驗準備認真閱讀實驗?zāi)康呐c實驗材料后要正確地解讀實驗，在此基礎(chǔ)上制定實驗計劃(修改、補充或編寫程序，提出實驗思

10、路，明確實驗步驟)，為上機實驗做好準備。四、實驗思路提示實驗步驟1、求解問題(2)中的食餌模型的微分方程組，并畫出解的圖形和相位圖。(1)以x=800, y=100為初始值，計算 x (t) ,y (t),當(dāng)t 0, 14時的數(shù)據(jù)。繪出解的圖形, 并分析捕食者和被捕食者的數(shù)量變化規(guī)律?？梢韵扔孟旅娴拿钋蠼馕鼋猓篊learx,y,txy=DSolvext=4*xt*xt*yt,yt=*yt+*xt*yt,x0=800,y0=100,x,y,t注：可以發(fā)現(xiàn)不能求出解析解。修改代碼如下，可以求數(shù)值解：Clearx,y,txy=NDSolvext=4*xt*xt*yt,yt=*yt+*xt*yt,x

11、0=800,y0=100,x,y,t,0,14繪出解的圖形：y=y/.xy1;x=x/.xy1;Plotxt,yt,t,0,14,PlotStyle-RGBColor0,0,1,RGBColor1,0,0圖捕食者和被捕食者的數(shù)量變化(2)以x為橫坐標，y為縱坐標繪制相位圖。根據(jù)圖形分析被捕食者數(shù)量增加 量的影響。繪制相位圖的命令：ParametricPlotxt,yt,t,0,14(減少)對捕食者數(shù)(M)圖相位圖2、用MATLAB求解問題(1)中Lorenz模型的微分方程。(1)打開MATLAB的編輯器；m文件中(2)在編輯器中用下面的幾個語句描述微分方程，并將其保存在的f unction x

12、dot = lorenzeq(t,x)xdot=k8/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);繪制出時間曲線與相空間曲線，如下圖所示。圖時間曲線與相空間曲線思考問題1、運用Mathematica求解Lorenz模型的微分方程組，從而了解系統(tǒng)狀態(tài)是如何變化的。2、求解以下問題(廣告的效用)：某公司生產(chǎn)一種耐用