微積分下余7.常微分方程

1、微分方程與計(jì)算機(jī)模擬 常微分方程數(shù)值解方法捕食者與被捕食者問(wèn)題有阻力拋射曲線問(wèn)題衛(wèi)星軌道模擬問(wèn)題數(shù)值方法求常微分方程初值問(wèn)題 求解步驟:(1)用函數(shù)文件定義一階微分方程(或方程組)右端函數(shù)；(2)用MATLAB命令ode23()求數(shù)值解或繪積分曲線。使用格式：T,Y = ode23(F,Tspan,y0)其中,Tspan = t0，tN是常微分方程求解區(qū)域，y0是初始值,F 是包括函數(shù)文件名字的符串。返回值(T,Y) 是求解區(qū)域內(nèi)離散數(shù)據(jù)及對(duì)應(yīng)數(shù)值解。例7.1 馬爾薩斯模型,以1994 年我國(guó)人口為12億為初值，求解常微分方程 N(t)表示人口數(shù)量,取人口變化率r =0.015,微分方程fun

2、ction z=fun1(t,N)z=0.015*N; ode23(fun1,1994,2020,12)T,N=ode23(fun1,1994,2020,12)命令窗口 編輯窗口 常微分方程組初值問(wèn)題一階常微分方程組初值問(wèn)題數(shù)值求解方法T,y = ode23( F ,Tspan,y0) 其中, F是函數(shù)文件, 表示 微分方程右端函數(shù)Tspan = t0 Tfinal 求解區(qū)域; y0 初始條件注: 函數(shù)F(t,y) 必須返回列向量. 數(shù)值解 y 的每一行對(duì)應(yīng)于列向量T中的每一行數(shù)據(jù)捕食者與被捕食者問(wèn)題 海島上有狐貍和野兔,當(dāng)野兔數(shù)量增多時(shí)，狐貍捕食野兔導(dǎo)致狐群數(shù)量增長(zhǎng);大量兔子被捕食使狐群進(jìn)入

3、饑餓狀態(tài)其數(shù)量下降;狐群數(shù)量下降導(dǎo)致兔子被捕食機(jī)會(huì)減少,兔群數(shù)量回升。微分方程模型如下計(jì)算 x(t)，y(t) 當(dāng)t0，20時(shí)的數(shù)據(jù)。繪圖并分析捕食者和被捕食者的數(shù)量變化規(guī)律。x(0)= 100y(0)=20 創(chuàng)建MATLAB的函數(shù)文件function z=fox(t,y)z(1,:)=y(1)-0.015*y(1).*y(2);z(2,:)=-y(2)+0.01*y(1).*y(2);Y0=100,20;t,Y=ode23(fox,0,20,Y0);x=Y(:,1);y=Y(:,2);figure(1),plot(t,x,b,t,y,r)figure(2),plot(x,y) 求微分方程數(shù)值

4、解并繪解函數(shù)圖形-兔子數(shù)量; -狐貍數(shù)量兔-狐數(shù)量變化相位圖拋射曲線實(shí)驗(yàn)，假設(shè)阻力與速度成正比。在微分方程中增加阻力項(xiàng) 符號(hào)計(jì)算方法syms t v g alfa kx=dsolve(D2x=-k*Dx,x(0)=0,Dx(0)=v*cos(alfa);y=dsolve(D2y=-g-k*Dy,y(0)=0,Dy(0)=v*sin(alfa);X=taylor(x,3,t),Y=simplify(taylor(y,3,t)2008電影集結(jié)號(hào)展現(xiàn)出視聽震撼的戰(zhàn)爭(zhēng)場(chǎng)面，92式山炮，炮彈初速:198米/秒,最大射程:2788米利用實(shí)驗(yàn)程序確定阻力系數(shù) kfunction Xmax=mlab72(k)

5、alfa=pi/4; v=198;g=9.8;t=0;dt=.1;x=0;y=0;while y=0; t=t+dt; xk=v*cos(alfa)*t-1/2*v*cos(alfa)*k*t2; yk=v*sin(alfa)*t-1/2*g*t2-1/2*t2*v*sin(alfa)*k; x=x,xk;y=y,yk;endXmax=xk;plot(x,y,ro)實(shí)驗(yàn)數(shù)據(jù):k 0.1 0.01 0.02 0.015Xmax 677.35 3073.15 2433.66 2719.33k=0.02k=0.015嫦娥一號(hào)軌道數(shù)據(jù)實(shí)驗(yàn)經(jīng)歷四次變軌提速后，衛(wèi)星才進(jìn)入地月轉(zhuǎn)移軌道。第一次變軌衛(wèi)星由初始

6、軌道進(jìn)入16小時(shí)軌道;第二次變軌衛(wèi)星進(jìn)入24小時(shí)軌道;第三次變軌衛(wèi)星進(jìn)入48小時(shí)軌道;第四次變軌衛(wèi)星進(jìn)入116小時(shí)地月轉(zhuǎn)移軌道。嫦娥一號(hào)衛(wèi)星進(jìn)入的初始軌道是周期為16小時(shí)的地球同步軌道。衛(wèi)星進(jìn)入初始軌道時(shí),最大速度大約為10.3(km/s),而奔月速度需要10.9(km/s)。假設(shè)五個(gè)軌道上最大速度從10.3(公里/秒)逐步增加到10.9(公里/秒)10.3，10.45，10.6，10.75，10.9根據(jù)牛頓萬(wàn)有引力定律，地球?qū)πl(wèi)星的引力大小為 地球引力參數(shù)：GM=3.986005105（km3/s2）衛(wèi)星運(yùn)動(dòng)方程轉(zhuǎn)換為一階微分方程組初始條件右端函數(shù)的函數(shù)文件function z=orbit(t,y)GM=3.986005e05;z(1,:)=y(2);z(2,:)=-GM*y(1)./(y(1).2+y(3).2).(3/2);z(3,:)=y(4);z(4,:)=-GM*y(3)./(y(1).2+y(3).2).(3/2);function Vmax,H=orbitlab(v,h,T)T0=T*60*60;Y0=-(6378+h),v*cos(-pi/2),0,v*sin(-pi/2);T,Y=ode23(orbit,0,T0,Y0);x=Y(:,1);y=Y(:,3);vx=Y(:,2);vy=Y(:,4);V=sqrt(v