CWG-微分方程求解Word版

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

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

3、方程 dsolve ('D2y+4*Dy+5*y=0','x')運行結(jié)果： ans =C1*exp(-2*x)*sin(x)+C2*exp(-2*x)*cos(x)即： 例3：解方程1 / 11 dsolve ('(1+x2)*Dy+2*x*y=x*exp(x2)','x')運行結(jié)果：ans =(1/2*exp(x2)+C1)/(1+x2)即： 如果要求微分方程的初值問題：，可輸入以下語句 dsolve ('D2y+4*Dy-2*y=0','y(0)=6','Dy(0)=10',&#

4、39;x')運行結(jié)果： ans = (3+11/6(1/2)*exp(-2+6(1/2)*x)+(-11/6(1/2)+3)*exp(-(2+6(1/2)*x)即： 2求微分方程數(shù)值解。求微分方程數(shù)值解命令為ode45，ode23，ode15s。對于不可以用積分方法求解的微分方程初值問題，可以用ode45，ode23，ode15s命令求特解。例4：求微分方程的近似解（）可用下面的命令： function f=odefun1(x,y) f=y2*x+y+x2; x,y=ode45(odefun1,0,4,3); plot ( x , y , r - )輸出結(jié)果為： 例5：解初值問題 ds

5、olve('Dy=-y+t+1','y(0)=1')輸出結(jié)果：ans =t+exp(-t)即： 現(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 off輸出結(jié)果：例6：求初值問題解：設(shè)，則原方程可化為Matlab語言：function f=odefun3(x,y)f=y(2);-y

6、(1)-sin(2*x);x,y=ode45('odefun3',pi,2*pi,1,1);plot(x,y(:,1),'r-')輸出結(jié)果：利用ode45命令還可以求解耦合微分方程，所謂耦合微分方程，方程組中的未知函數(shù)是相互影響的，相互依賴的，其中的一個求解會影響到另一個求解，下面求一對耦合微分方程的數(shù)值解：例7：解方程組其中Matlab語句： function f=odefun4(t,y)f=y(2),-0.01*y(2)-sin(y(1)'t,y=ode15s('odefun4',0,100,0,2.1);函數(shù)的圖像：plot(t,y

7、(:,1),'r-')輸出結(jié)果為：函數(shù)的圖像：plot(t,y(:,2),'r-')輸出結(jié)果：用生成參數(shù)圖形plot(y(:,1),y(:,2),'r-')輸出結(jié)果：五、 實驗過程1歐拉折線法對于初值問題，我們考慮函數(shù)的線性近似 由于函數(shù)可微，在包含的一個很小的鄰域內(nèi)是得很好的近似。歐拉折線法就是通過一系列的線性近似得到在較大區(qū)間內(nèi)的的近似解。第一步：設(shè)，其中很小，則 是得很好的近似，在區(qū)間（無妨設(shè)）上能被很好的近似。第二步：利用和斜率來進行下一步近似，設(shè)，由 近似表示。第三步：利用點和斜率，對于，由 近似表示。 這樣我們就得到一列點列。而連接這

8、個點列的折線就是初值問題的一個近似解。這就是所謂的歐拉折線法。其一般的步驟是 例8：利用歐拉折線法球初值問題的近似解。以下是求此初值問題的Matlab語句 function odefun6(n,d) X=0,1; for k=1:n/d X(k+1,1)=X(k,1)+d; X(k+1,2)=X(k,2)+(1+X(k,2)*d; end plot(X(:,1),X(:,2) odefun6(1,0.01)此初值問題的精確解為，以上語句可以實現(xiàn)對精確解和近似解的圖像進行比較。 hold on x=0:0.01:1; y=2*exp(x)-1; plot(x,y,'r-')hol

9、d off輸出結(jié)果：2微分方程的斜率場例9：求一階微分方的斜率場。一階微分方程求斜率場的Matlab語句如下： function odefun7(inx,axx,iny,axy) a=(axx-inx)/0.1; b=(axy-iny)/0.1; z=inx,iny;m=1;s=1; for j=1:b for 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; end s=m; end for i=1:length(z) x=z(i,1);y=z(i,2); x1=x+0.1;y1=y+(y2*x+y)*0.1; plot(x

10、,x1,y,y1,'b-') plot(z(:,1),z(:,2),'b.','Markersize',2) hold on end hold off odefun7(-2,2,-2,2)輸出結(jié)果為六、 結(jié)論與應(yīng)用研究衛(wèi)星繞地球運行的軌跡。根據(jù)牛頓第二運動定律：和萬有引力定理。所以，其中為地球的質(zhì)量，為衛(wèi)星所在位置的坐標(biāo)，。因此我們有在軸上加速度分量為，在軸上加速度分量為，設(shè)衛(wèi)星的運動方程為，則有。如果我們假定衛(wèi)星以初速度在處入軌，地球質(zhì)量為。，這是一個初值問題。設(shè)微分方程可化為：Matlab語句：function f=odefun5(t,y,f

11、lag,G,M)r=sqrt(y(1)2+y(2)2);f=y(3),y(4),-G*(M/r3)*y(1),-G*(M/r3)*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-')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','none');camlight right; lighting phonghold off輸出結(jié)果：七、 練習(xí)1 求下列微分方程的通解。1）2）3）2 求初值問題的解。3 求微分方程在初始條件下的精確解和用歐拉折線法得到的近似解并作圖。注意觀察的選取對解的精確度的影響。4 求微分方程組的特解，并作出用生成參數(shù)圖形。5 請用歐拉折線法在Matla