matlab：常微分方程的應(yīng)用

實(shí)驗(yàn)名稱：常微分方程的應(yīng)用姓名：學(xué)號(hào)： 實(shí)驗(yàn)日期：年5月10日MATLAB有關(guān)求解常微分方程（組）驗(yàn)介紹常微分方程（組）數(shù)值求解Euler彈系統(tǒng)改進(jìn)和產(chǎn)品銷售量等幾個(gè)實(shí)際問題的微分方程建模及求解方法。實(shí)驗(yàn)項(xiàng)目：1、求解微分方程組，并繪制Appolo衛(wèi)星軌跡圖2、廣告效應(yīng)實(shí)驗(yàn)背景：1、除常系數(shù)線性微分方程可用特征根法求解，少數(shù)特殊方程可用初等積分法求解外，大部分微分方程無顯示解，應(yīng)用中主要依靠數(shù)值解法。2、在商業(yè)活動(dòng)中，常常會(huì)需要制定銷售策略來預(yù)測(cè)公司產(chǎn)品的銷售量。這就需要建立數(shù)學(xué)模型來預(yù)測(cè)銷量。實(shí)驗(yàn)具體過程：1、題目：已知Appolo衛(wèi)星的運(yùn)動(dòng)軌跡（x，y）滿足下面的方程d2x dy x）x）=2 +x- - ，dt2

dt 3 31 2d2y=-2dx+y-y-y，dt2 dt 3 31 2其中=1/82.45,=1-1

=x2+y22

=x2+y2（0=1.，=0（0）=’004935371Appolo衛(wèi)星軌跡圖。理論推導(dǎo)或編程說明：首先引入輔助變量t=t，y1

=x，y2

=dx，ydt 3

=y，y4

=dy，dt化為一階方程組dy1=ydt 2

dy

2=2y+y

（y1

）（y-）- 1dt 4 1dy3=ydt 4

3 31 2dy y y 4=-2ydt

+y3

3-3331 233先寫M函數(shù)p113_4.mfunctionf=p113_4fun（t，y）

%M函數(shù)p113_4.mf=[ y(2) ; 2*y(4)+y(1)-（+）/

3（）/1

3 ; y(4) ;2-2*y(2)+y(3)-*y(3)/

3-*y(3)/3]1 2=1/82.45,=1-,1再在指令窗口執(zhí)行

=2

，=2+22>>y0=[0,0,0，-1.04935371]；>>[t，y]=ode4（@p114_4fu，[024]，y；>>plo（（，1，（，3）；程序（或命令：M函數(shù)p113_4fun.mM函數(shù)p113_4fun.mfunctionf=p113_4fun(t,y)u=1/82.45;l=1-u;r(1)=sqrt((y(1)+u)^2+y(3)^2);r(2)=sqrt((y(1)+l)^2+y(3)^2);f=[y(2);2*y(4)+y(1)-l*(y(1)+u)/r(1)^3-u*(y(1)-l)/r(2)^3;y(4);-2*y(2)+y(3)-l*y(3)/r(1)^3-u*y(3)/r(2)^3];再在指令窗口執(zhí)行>>y0=[1.2,0,0,-1.04935371];>>[t,y]=ode45(@p113_4fun,[024],y0);>>plot(y(:,1),y(:,3));對(duì)實(shí)驗(yàn)題目的解答：改進(jìn)程序或思考：clear;dfun=inline('[x(2);2*x(3)+x(1)-((1-1/82.45)*(x(1)+1/82.45))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*(x(1)-1+1/82.45))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3;x(4);-2*x(2)+x(3)-((1-1/82.45)*x(3))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*x(3))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3]','t','x');[t,x]=ode45(dfun,[024],[1.2;0;0;-1.04935371]);plot(x(:,1),x(:,3));2、題目：(廣告效應(yīng))某公司生產(chǎn)一種耐用消費(fèi)品，市場(chǎng)占有率為5%時(shí)開始做廣0.5。建立該問題的數(shù)學(xué)模型，分別求其解析解和數(shù)值解，并作比較；理論推導(dǎo)或編程說明：x'由模型x

0.5(1x)

dxdt

0.5(xx2)2 dxdtx)2(1 1 )dxdt兩邊同時(shí)積分，得

x 1x2x

(1 1 )dxtdt0.05

x 1x 0t2(lnxx)ln0.05ln序（或命令：>>fun=inline('0.5*(1-x)*x','t','x');>>[t,x]=ode45(fun,[010],0.05)；>>[t,x]ans=00.05000.10580.05260.21150.05530.31730.05810.42310.06110.67310.06860.92310.07711.17310.08641.42310.09681.67310.10831.92310.12102.17310.13492.42310.15022.67310.16692.92310.18503.17310.20463.42310.22573.67310.24833.92310.27234.17310.29784.42310.32464.67310.35254.92310.38165.17310.41155.42310.44205.67310.47315.92310.50436.17310.53556.42310.56646.67310.59686.92310.62657.17310.65527.42310.68297.67310.70937.92310.73448.17310.75818.42310.78028.67310.80098.92310.82019.17310.83789.42310.85419.56730.86299.71150.87129.85580.8790>>plot(x,t,'.')；>>holdon;>>fplot(inline('2*(log(x)-log(1-x)-log(0.05)+log(0.95))'),[0.05,0.9])>>holdoff;>>id=min(find(x>0.8));>>t(id)ans=8.6731對(duì)實(shí)驗(yàn)題目的解答：解：x' x'xx

0.5(1x)，數(shù)值解可用MATLAB求解：>>fun=inline('0.5*(1-x)*x','t','x');>>[t,x]=ode45(fun,[010],0.05)；>>[t,x]ans=00.05000.10580.05260.21150.05530.31730.05810.42310.06110.67310.06860.92310.07711.17310.08641.42310.09681.67310.10831.92310.12102.17310.13492.42310.15022.67310.16692.92310.18503.17310.20463.42310.22573.67310.24833.92310.27234.42310.32464.67310.35254.92310.38165.17310.41155.42310.44205.67310.47315.92310.50436.17310.53556.42310.56646.67310.59686.92310.62657.17310.65527.42310.68297.67310.70937.92310.73448.17310.75818.42310.78028.67310.80098.92310.82019.17310.83789.42310.85419.56730.86299.7115