CWG-微分方程求解

1、微分方程求解實(shí)驗(yàn)?zāi)康呐c要求1掌握用 Matlab 求微分方程及其方程組解的方法； 2學(xué)會(huì)求微分方程近似解的歐拉折線法； 3學(xué)會(huì)建立一些簡(jiǎn)單問(wèn)題的微分方程模型，并能運(yùn)用Matlab 分析研究這些問(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í)，我們一般須先建立微分方程，再利

2、用 所學(xué)的數(shù)學(xué)知識(shí)解微分方程。 事實(shí)上真正能找到精確解的微分方程只是很少一部 分，大局部只能求近似解，即數(shù)值解。四、背景知識(shí)介紹1 求微分方程解析解的命令。 求微分方程解析解的命令為：dsolve'方程1''方程2'，初始條件1''初始條件2'，自變量 對(duì)于可用積分方法求解的微分方程和微分方程組，可以用dsolve命令來(lái)求其通解和特解。例 1：要求方程 y 3y 4y 0的通解，可以輸入以下語(yǔ)句 Matlab 命令： dsolve (D2y+3*Dy-4*y=0 ','x ')運(yùn)行結(jié)果：ans =C1*exp(-4

3、*x)+C2*exp(x)D2 表示，三階導(dǎo)數(shù)用 't'。D3 表示，以此類yC1e 4xC2ex注：求一階用 D 表示，二階導(dǎo)數(shù)用 推。如果自變量沒(méi)有選定，默認(rèn)自變量為例 2：解方程 y4y 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)即：y e 2x C1 sin x C2 cosx例 3：解方程 1 x2 y 2xy xex2dsolve ('(1+xA2)*Dy+2*x*y=x*exp(xA2)',

4、9;x')即:運(yùn)行結(jié)果：ans =(1/2*exp(xA2)+C1)/(1+xA2) lex2 G y V1 x如果要求微分方程的初值問(wèn)題：y 4y 2y 0x 06 , y|x 010，可輸入以下語(yǔ)句dsolve ('D2y+4*Dy-2*y=0','y(0)=6','Dy(0)=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,

5、 ode15so對(duì)于不可以用積分方法求解的微分方程初值問(wèn)題，可以用ode45, ode23,ode15s命令求特解。例4:求微分方程yy2x y x2 , y x 0 3的近似解0 x 4可用下面的命令：fun cti onf=odef un 1(x,y)f=yA2*x+y+xA2;x,y=ode45( 6defun1 '0,4,3);plot ( x , y , r -'輸出結(jié)果為：例5:解初值問(wèn)題y y t 1, y 01dsolve('Dy=-y+t+1','y(0)=1')輸出結(jié)果：ans =t+exp(-t)即：y t e t現(xiàn)在我們用

6、數(shù)值求解命令求解后和解析解比擬 function f=odef un 2(t,y)f=-y+t+1;t=0:0.1:1; y=t+exp(-t); plot(t,y,'b-') hold ont,y=ode45('odefu n2',0,1,1); plot(t,y,'r.')hold off輸出結(jié)果：例6:求初值問(wèn)題y y sin 2x 0 , y 1 , y 1解：設(shè)y1y , y2y，那么原方程可化為y1 y2y2y1 sin 2xy11 , y21Matlab 語(yǔ)言：fun cti on f=odef un 3(x,y)f=y(2);-y

7、(1)-si n(2*x);x,y=ode45('odefu n3',pi,2*pi,1,1);plot(x,y(:,1),'r-')輸出結(jié)果:利用ode45命令還可以求解耦合微分方程，所謂耦合微分方程，方程組中 的未知函數(shù)是相互影響的，相互依賴的，其中的一個(gè)求解會(huì)影響到另一個(gè)求解， 下面求一對(duì)耦合微分方程的數(shù)值解：例7:解方程組xt yt其中x 00 , y 02.1y t 0.01y t sin x tMatlab 語(yǔ)句：function f=odef un 4(t,y) f=y (2),-0.01*y (2)-si n( y(1)' t,y=ode

8、15s('odefu n4',0,100,0,2.1);函數(shù)x x t的圖像：plot(t,y(:,1),'r-')輸出結(jié)果為：函數(shù)y y t的圖像:plot(t,y(:,2),'r-')輸出結(jié)果:用xt , yt生成參數(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可微，在包含X0的一個(gè)很小的鄰域內(nèi)L x是y x得很好的近第一步 ：設(shè) X1 X0X ，其中 X 很小，

9、那么y1 L X1y0y X0X1X0y0 f X0 , y0 X1X0是y Xi得很好的近似，在區(qū)間X0 ,X1 無(wú)妨設(shè) X 0 上 y X能被 L X 很好的似。歐拉折線法就是通過(guò)一系列的線性近似得到在較大區(qū)間內(nèi)的y x 的近似解近似。第二步 ：利用 x1 , y1 和斜率 f x1 , y1 來(lái)進(jìn)行下一步近似，設(shè) x2 x1 x ， y x2 由y2 y1 f x1 , y1 x2 x1近似表示。第三步 ：利用點(diǎn) x2 , y2 和斜率 f x2 , y2 ，對(duì)于 x3 x2 x ， y x3 由 y3 y2 f x2 , y2 x3 x3近似表示。這樣我們就得到一列點(diǎn)列 x0 , y0

10、 , x1 , y1 , x2 , y2 , x3 , y3 。而連 接這個(gè)點(diǎn)列的折線就是初值問(wèn)題 y f x, y , y x0 y0 的一個(gè)近似解。這就是 所謂的歐拉折線法。其一般的步驟是X1X0dX ,y1y0f X0 , y0 X1X0X2X1dX ,y2y1f X1 , y1 X2X1X3X2dX ,y3y2f X2 , y2 X3X2XnXn1 dX, ynyn1f Xn 1 , yn1 XnXn 1例 8：利用歐拉折線法球初值問(wèn)題 y 1 y , y 0 1 的近似解 以下是求此初值問(wèn)題的 Matlab 語(yǔ)句function odefun6(n,d)X=0,1;for k=1:n

11、/dX(k+1,1)=X(k,1)+d;X(k+1,2)=X (k,2)+(1+X(k,2)*d;endplot(X(:,1),X(:,2)odefu n6(1,0.01)此初值問(wèn)題的精確解為y2ex 1，以上語(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：求一階微分方y(tǒng)y2x y的斜率場(chǎng)。一階微分方程求斜率場(chǎng)的Matlab語(yǔ)句如下：fun ctio n odef un 7(i nx,axx,i ny ,axy)a=(axx-i nx)/0

12、.1;b=(axy-i ny)/0.1;z=i nx,i ny ;m=1;s=1;for j=1:bfor k=1:am=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:le ngth(z)x=z(i,1);y=z(i,2);x1=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 offodefu n7

13、(-2,2,-2,2)輸出結(jié)果為30卜-21£500.511.622.6六、結(jié)論與應(yīng)用研究衛(wèi)星繞地球運(yùn)行的軌跡根據(jù)牛頓第二運(yùn)動(dòng)定律:F mad2rm2dt和萬(wàn)有引力定理FmMG-。所以ra G耳，其中M為地球的質(zhì)量，x, yr因此我們有在x軸上加速度分量為a為衛(wèi)星所在位置的坐標(biāo)，r x2 y2。Gx，r在y軸上加速度分量為ayGMry，設(shè)衛(wèi)星的運(yùn)動(dòng)方程為Xrxtyt，那么有MG飛xr。如果我們M gF r假定衛(wèi)星以初速度vy 04000m/s 在 x 04.2 107m處入軌，地球質(zhì)量為5.97 1024kg。G 6.672 10 11 N m2/kg2 ，這是一個(gè)初值問(wèn)題。設(shè)x ,

14、 y4 y微分方程可化為:y2y4y3y4gMryigMry2Matlab 語(yǔ)句:yi 04.2 107 , y2 00, yi 00, y? 04000function f=odefun5(t,y,flag,G ,M) r=sqrt(y(1)A2+y (2) A2); 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,G ,M); plot(y(:,1),y(:,2),'r-&#

15、39;)hold onX,Y,Z=sphere(10);axis('image')R=0.64e7;X=R*X;Y=R*Y;z=0*Z;surf(X,Y,Z,'FaceColor','red','EdgeColor',' non e');camlightright;phonghold off輸出結(jié)果：lighti ng七、練習(xí)1 求以下微分方程的通解。1y 42yy02y2y5yex sin 2x3y6y9y彳 3xx 1 e2.求初值問(wèn)題y ysin2x 0 , y x 1 , y x 1 的解。3 求微分方程x21y 2xy cosx 0在初始條件y%。1下的精確解和用歐拉折線法得到的近似解0 x 1并作圖