CWG-微分方程求解

1、微分方程求解實(shí)驗(yàn)?zāi)康呐c要求1掌握用Matlab求微分方程及其方程組解的方法；2學(xué)會(huì)求微分方程近似解的歐拉折線法；Matlab分析研究這些3學(xué)會(huì)建立一些簡(jiǎn)單問(wèn)題的微分方程模型，并能運(yùn)用 問(wèn)題。二、問(wèn)題描述對(duì)于很多實(shí)際問(wèn)題，要直接找出所需的函數(shù)關(guān)系往往非常困難，但根據(jù)實(shí) 際問(wèn)題所提供的條件，有時(shí)卻可以列出含有未知函數(shù)導(dǎo)數(shù)的關(guān)系式， 這樣的關(guān)系 式就是所謂的微分方程。怎樣利用微分方程求得所需未知函數(shù)， 往往是我們解決 實(shí)際問(wèn)題經(jīng)常需要面對(duì)的問(wèn)題，即解微分方程。這里我們借用Matlab對(duì)此問(wèn)題進(jìn)行簡(jiǎn)單探討。三、問(wèn)題分析在處理關(guān)于微分方程的實(shí)際問(wèn)題時(shí)，我們一般須先建立微分方程，再利用 所學(xué)的數(shù)學(xué)知識(shí)解微

2、分方程。事實(shí)上真正能找到精確解的微分方程只是很少一部 分，大部分只能求近似解，即數(shù)值解。四、背景知識(shí)介紹1.求微分方程解析解的命令。求微分方程解析解的命令為：dsolve（方程1方程2，初始條件1初始條件2，自變量），對(duì)于可用積分方法求解的微分方程和微分方程組，可以用dsolve命令來(lái)求其通解和特解。例1:要求方程y”的通解，可以輸入以下語(yǔ)句 Matlab命令:運(yùn)行結(jié)果:即dsolve (D2y+3*Dy-4*y=0 ,) ans =C1*exp(-4*x)+C2*exp(x)4xxy 二 CeC?e注：求一階用D表示，二階導(dǎo)數(shù)用D2表示，三階導(dǎo)數(shù)用D3表示，以此類(lèi) 推。如果自變量沒(méi)有選定，默

3、認(rèn)自變量為 ）。例2:解方程y “ 4y * 5y = 0dsolve (D2y+4*Dy+5*y=0,x)運(yùn)行結(jié)果：ans =C1*exp(-2*x)*sin(x)+C2*exp(-2*x)*cos(x)例 3:解方程 1 x2 y 2xy = xedsolve (1+xA2)*Dy+2*x*y=x*exp(xA2),x)運(yùn)行結(jié)果：ans =(1/2*exp(xA2)+C1)/(1+xA2)1 2-exG即：y = 2y1 x2如果要求微分方程的初值問(wèn)題：y4y_2y=0 y x = 6 , yx衛(wèi)=1。，可輸入以下語(yǔ)句dsolve (D2y+4*Dy-2*y=0,y(0)=6,Dy(0)=

4、10,x)運(yùn)行結(jié)果：ans =(3+11/6A(1/2)*exp(-2+6A(1/2)*x)+(-11/6A(1/2)+3)*exp(-(2+6A(1/2)*x)即:2 求微分方程數(shù)值解。求微分方程數(shù)值解命令為 ode45, ode23, ode15s對(duì)于不可以用積分方法求解的微分方程初值問(wèn)題，可以用ode45, ode23ode15s命令求特解。例4：求微分方程y = y2x十y+x2 , y x=0 = 3的近似解（0Ex蘭4）可用下面的命令:function f=odefun1(x,y) f=yA2*x+y+xA2; x,y=ode45( 6defun1,0,4,3);plot ( x

5、, y , r -)輸出結(jié)果為:3.7例5:解初值問(wèn)題y - -y 71, yO=1dsolve(Dy=-y+t+1,y(0)=1) 輸出結(jié)果：ans =t+exp(-t)即：y = t e丄現(xiàn)在我們用數(shù)值求解命令求解后和解析解比較function f=odefun2(t,y)f=-y+t+1;t=0:0.1:1; y=t+exp(-t); plot(t,y,b-) hold ont,y=ode45(odefun2,0,1,1); plot(t,y,r.) hold off1 4例 6:求初值問(wèn)題 y y sin 2x =0 , yi： 1 , yi： 1解：設(shè)yy ,，則原方程可化為yi =

6、 y2* y2 = _yi _sin 2xyS )=i , y )=iMatlab 語(yǔ)言：function f=odefun3(x,y) f=y(2);-y(1)-sin(2*x); x,y=ode45(odefun3,pi,2*pi,1,1); plot(x,y(:,1),r-)1 544.5G65利用ode45命令還可以求解耦合微分方程，所謂耦合微分方程，方程組中 的未知函數(shù)是相互影響的，相互依賴的，其中的一個(gè)求解會(huì)影響到另一個(gè)求解， 下面求一對(duì)耦合微分方程的數(shù)值解： 例7：解方程組恍Txt )其中心。，心1Matlab 語(yǔ)句：function f=odefun4(t,y) f=y(2),

7、-0.01*y(2)-sin(y(1); t,y=ode15s(odefun4,0,100,0,2.1);函數(shù)x =x t的圖像：plot(t,y(:,1),r-)輸出結(jié)果為:26函數(shù)y = y t的圖像：plot(t,y(:,2),r-)輸出結(jié)果：用xtl, yt 1生成參數(shù)圖形plot(y(:,1),y(:,2),r-)輸出結(jié)果:五、實(shí)驗(yàn)過(guò)程1 歐拉折線法對(duì)于初值問(wèn)題y丄f x, y , y X。二y，我們考慮函數(shù)y x的線性近似Lx 二 y x。 y x。x - x。由于函數(shù)y x可微，在包含X。的一個(gè)很小的鄰域內(nèi)L x是y x得很好的近 似。歐拉折線法就是通過(guò)一系列的線性近似得到在較大

8、區(qū)間內(nèi)的yx的近似解。第一步：設(shè)x1 =x。r：x，其中x很小，貝Uyi = LXi二y。y x。Xi - x。二y。fx。，y。Xi - x。是y Xi得很好的近似，在區(qū)間X。，Xi丨（無(wú)妨設(shè)x 0） 上 y x能被L x很好的 近似。第二步：利用x-i , yi和斜率f治，yj來(lái)進(jìn)行下一步近似，設(shè)xx x， y X2由y2 二 yi f Xi , yi X2 - Xi近似表示。第三步：利用點(diǎn)x2 , y2和斜率f x2 , y2，對(duì)于X3 = x2 ： =x，y x3由yY2 f X2 , y2 X3 -X3近似表示。這樣我們就得到一列點(diǎn)列Xo ,yo , Xi,yi ,X2, y2 ,

9、X3 , y3。而連接這個(gè)點(diǎn)列的折線就是初值問(wèn)題yfx,y,y Xo=y0的一個(gè)近似解。這就是所謂的歐拉折線法。其一般的步驟是XidX ,yi=y fX。，yXi - XoX2 = xi dx ,y2= yifxi ,yiX2 - XiX3 = X2 dx ,y3= y2fX2 ,y2X3 - xXn = Xnj dX , yn = yn f X.,y.X. - X.例8：利用歐拉折線法球初值問(wèn)題 yi y , yO =i的近似解。以下是求此初值問(wèn)題的Matlab語(yǔ)句function odefun6(n,d)X=0,1;for k=1:n/dX(k+1,1)=X(k,1)+d;X(k+1,2)

10、=X(k,2)+(1+X(k,2)*d;endplot(X(:,1),X(:,2)odefun6(1,0.01)此初值問(wèn)題的精確解為y=2ex-i，以上語(yǔ)句可以實(shí)現(xiàn)對(duì)精確解和近似解的圖像進(jìn)行比較。hold onx=0:0.01:1;y=2*exp(x)-1;plot(x,y,r-)hold off輸出結(jié)果：2 微分方程的斜率場(chǎng)例9：求一階微分方目二y2x y的斜率場(chǎng)。一階微分方程求斜率場(chǎng)的Matlab語(yǔ)句如下：function odefun7(inx,axx,iny,axy) a=(axx-inx)/0.1;b=(axy-iny)/0.1;z=i nx,i ny;m=1;s=1;for j=1

11、:bfor k=1:a m=m+1;z(k+s,1)=z(k,1)+0.1;z(k+s,2)=z(k,2)+(j-1)*0.1; ends=m;endfor i=1:length(z)x=z(i,1);y=z(i,2);x仁 x+0.1;y 1=y+(yA2*x+y)*0.1; plot(x,x1,y,y1,b-) plot(z(:,1),z(:,2),b.,Markersize,2) hold onendhold offodefun7(-2,2,-2,2)輸出結(jié)果為六、結(jié)論與應(yīng)用研究衛(wèi)星繞地球運(yùn)行的軌跡根據(jù)牛頓第二運(yùn)動(dòng)定律：2F -ma - m 2和萬(wàn)有引力疋理F - -G 2。所以dtra

12、二-，其中M為地球的質(zhì)量，x , y為衛(wèi)星所在位置的坐標(biāo)，r二x2 y2 。 r因此我們有在x軸上加速 度分量為ax = -GM3x，在y軸上加速度分量為 ray二-G My，設(shè)衛(wèi)星的運(yùn)動(dòng)方程為rx= x(t),y = y(t)怙k c Mx (t )= -G 飛 x，則有rMy t = -G 飛 y Lr如果我們假定衛(wèi)星以初速度 vy 0 =4000m/s在x 0 - -4.2 107m處入軌，地球質(zhì)量為M -5.97 1024kg。G =6.672 101N m2/kg2 ，這是一個(gè)初值 問(wèn)題。設(shè)yi = x , y2 = y , y3 = x , y4 = y 微分方程可化為:y； =

13、y3y； =y4. mr-My4 = -G py2ry O A-4.2 107 , y； 01=0, y； 0A0, y； 01=4000Matlab 語(yǔ)句：function f=odefun5(t,y,flag,GM) r=sqrt(y八2+y(2)八2); f=y(3),y(4),-G*(M/rA3)*y(1),-G*(M/rA3)*y(2);G=6.672e-11;M=5.97e24; t,y=ode45(odefun5,0,60*60*24*6,-4.2e7,0,0,4000,GM);plot(y(:,1),y(:,2),r-)hold onX,Y,Z=sphere(10);axis(image)R=0.64e7;X=R*X;Y=R*Y;z=0* Z;lightingsurf(X,Y,Z,FaceColor,red,EdgeColor,none);camlight right; phong hold off輸出結(jié)果：七、練習(xí)1 .求下列微分方程的通解1) y 4 -2y y =02) y -2y 5y = ex sin 2x3) y 辭6y 9y = x 1 e3x2.求初值問(wèn)題 y + y + sin2x=0 , y =1 , y =1 的解3.求