1、學生實驗報告（4）學生姓名陳宋學號008888同組人：無實驗項目Matlab程序設計與作圖必修 選修演示性實驗驗證性實驗 操作性實驗 綜合性實驗實驗地點H110實驗儀器臺號D02指導教師謝建實驗日期及節次2011-04-18第八九十節、實驗綜述.歸納和學習求解常微分方程（組）的基本原理和方法;2. 掌握解析、數值解法，并學會用圖形觀察解的形態和進行解的定性分析;3. 熟悉MATLAB軟件關于微分方程求解的各種命令;4. 通過范例學習建立微分方程方面的數學模型以及求解全過程;通過該實驗的學習，使學生掌握微分方程（組）求解方法（解析法、歐拉法、梯 度法、改進歐拉法等），對常微分方程的數值解法有一個

2、初步了解，同時學會使用MATLAB軟件求解微分方程的基本命令，學會建立微分方程方面的數學模型。這對 于學生深入理解微分、積分的數學概念，掌握數學的分析思維方法，熟悉處理大量 的工程計算問題的方法是十分必要的。二、實驗過程（實驗步驟、記錄、數據、分析）1.開啟MATLAB軟件平臺，開啟 MATLAB編輯窗口;2根據問題，建立的線性規劃模型，并編寫求解規劃模型的 M文件;3.保存文件并運行；4 .觀察運行結果（數值或圖形），并不斷地改變參數設置觀察運行結果;5. 根據觀察到的結果和體會，寫出實驗報告。三、實驗要求與任務根據實驗內容和步驟，完成以下實驗，要求寫出實驗報告（實驗目的-問題-數學模型-算

3、法與編程-計算結果-分析、檢驗和結論)1.求微分方程的解析解，并畫出它們的圖形。y '= y + 2 x, y (0) = 1,0< x <1 ；程序如下：由 y=dsolve('Dy=y+2*x','y(0)=1','x')得出解析解 y =-2*x-2+3*exp(x)建立函數 m 文件：function y=myfun4(x)y=-2*x-2+3*ex p(x)畫圖函數為fplot('myfu n4',0,1)圖形如下:2.求微分方程首先建立函數7(0) = 0; u(0) = 0的數值解，要求編寫求解程

4、序。t =0 10M 文件：function dy=myfun5(t,y)dy=zeros(2,1);dy(i)=y(2)；dy( 2)=0.1*(y(1)A3)-y(1);輸入命令：T,Y=ode15s('myfu n5',0,10,0,1);plot(T,丫(:,1),'-',T,丫(:,2),'*');圖形如下:3. Rossler微分方程組:X = -y - ZI* y = X + ayIz = b 中 z(x - c)當固定參數b=2, c=4時，試討論隨參數a由小到大變化(如 a (0, 0.65)而方程解的變化情況，并且畫出空間曲線

5、圖形，觀察空間曲線是否形成混沌狀?首先建立函數M文件rossler.m ，在其中用x(1)表示x，用x(2)表示y，用x(3)表示乙function r=rossler(t,x) global a;global b;global c;r=-x (2)-x(3);x(1)+a*x(2);b+x(3)*(x(1)-c);主程序如下:global a;global b;global c;b=2;c=4;t0=0,200;for a=0:0.02:0.65t,x=ode45('rossler',t0,0,0,0);asub plot(1,2,1);plot(t,x(:,1),'

6、r',t,x(:,2),'g',t,x(:,3),'b');title('x(紅色),y(綠色),z(籃色)隨t變化情況');xlabel('t'); sub plot(1,2,2);plot3(x(:,1),x(:,2),x(:,3)title('相圖');xlabel('x');ylabel('y');zlabel('z');p ause end08,當 a=0 時，(x,y,z)收斂于(0,0.5,0.5) 當a=0.05時，(x,y,z)仍然收斂，但收

7、斂速度較小。其收斂速度已當a=0.12時，（x,y,z）仍然收斂（準確的說，只能說具有收斂的趨勢,大大降低）。法曲 為舟俸色）匿邑1砒舷04'0.2Q-0.3-06-060.6八 2(50100150200a=0.l20丄0 05a=0.16時，(x,y,z)已經發散。但是,(x,y,z)并不是發散于無窮大，而是周期性變化。40a=0.16220501001 麵 200411.a=U.2O15,05*12£501001507 -5a-0348O365 我也/沌）丄潘a刪變化AS說I.& Ij.Q50 IDO la 2004.炮彈發射角的確定炮彈發射視為斜拋運動，已知初

8、始速度為200m/s,若要擊中水平距離360m、垂直距離160m的目標，當忽略空氣阻力時，發射角應為多大？此時炮彈的運行軌 跡如何?要求:建立在忽略空氣阻力情況下的描述炮彈發射軌跡的數學模型;用Matlab軟件求解方程和微分方程;結合實際對解的合理性進行分析。進一步思考:如果要考慮水平方向的阻力，且設阻力與（水平方向）速度成正比，系數為0.1（ 1/s），結果又如何？此時炮彈的運行軌跡如何?建立模型x=200cos 0 *t2y=200sin0 *t-1/2*9.8t2函數myfun6.m文件如下:fun cti on y=myfun 6(x)y=360*ta n(x)-4.9*(360/20

9、0/cos(x)A2-160;建立Newton2.m函數如下:fun cti on y=n ewt on 2(x0,x1, n,tol)x(1)=x0;x(2)=x1;b=1;i=2;while(abs(b)>e ps*x(i)x(i+1)=x(i)-myfu n6(x(i)*(x(i)-x(i-1)/(myfu n6(x(i)-myfu n6(x(i-1);b=x(i+1)-x(i);i=i+1;if(i>n) error(' n is full');endenddis p(i-2);y=x(i);運行結果如下：>> nevton2(i0, 5, Ij

10、 lOOj le-5)7ans =0,4633建立模型如下:A噲代入初始條件可以得出 x=-10*200cos0*ex p(-0.1t)+10*200cos 0建立 myfun7 函數如下：function y=myfun7(x) y=200*sin(x)*(-10*log(1-360/2000/cos(x)-4.9*(-10*log(1-360/2000/cos(x).八2)-160建立newton2函數如下： fun cti on y=n ewt on 2(x0,x1, n,tol) x(1)=x0;x(2)=x1;b=1;i=2；while(abs(b)>e ps*x(i)x(i+1)=x(i)-myfu n7(x(i)*(x(i)-x(i-1)/(myfu n7(x(i)-myfu n7(x(i-1);b=x(i+1)-x(i);i=i+1;if(i> n) error('n is full'); endenddis p(i-2);y=x(i);運行結果如下：>> newton2(0. 5, 1, 100, le-6)441.4577441.457