MATLAB求解微分方程實(shí)驗(yàn)6

實(shí)驗(yàn)ExperimentsinMathematics微分方程求解實(shí)驗(yàn)?zāi)康膶?shí)驗(yàn)內(nèi)容MATLAB2、學(xué)會(huì)用Matlab求微分方程的數(shù)值解.實(shí)驗(yàn)軟件1、學(xué)會(huì)用Matlab求簡(jiǎn)單微分方程的解析解.1、求簡(jiǎn)單微分方程的解析解.2、求微分方程的數(shù)值解.微分方程的解析解求微分方程（組）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始條件’,‘自變量’)注意：①y'Dy，y''D2y②自變量名可以省略，默認(rèn)變量名‘t’。例1輸入：y=dsolve('Dy=1+y^2')y1=dsolve('Dy=1+y^2','y(0)=1','x')輸出：y=tan(t-C1)（通解）

y1=tan(x+1/4*pi)（特解）MATLAB軟件求解例2常系數(shù)的二階微分方程y=dsolve('D2y-2*Dy-3*y=0','x')y=dsolve('D2y-2*Dy-3*y=0','y(0)=1,Dy(0)=0','x')輸入：y=C1*exp(-x)+C2*exp(3*x)y=3/4*exp(-x)+1/4*exp(3*x)結(jié)果：x=dsolve('D2x-(1-x^2)*Dx+x=0','x(0)=3,Dx(0)=0')例3非常系數(shù)的二階微分方程無解析表達(dá)式！x=dsolve('(Dx)^2+x^2=1','x(0)=0')例4非線性微分方程x=sin(t)-sin(t)若欲求解的某個(gè)數(shù)值解，如何求解？t=pi/2;eval(x)MATLAB軟件求解輸入：[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y')[x,y]=dsolve('Dx=3*x+4*y','Dy=-4*x+3*y','x(0)=0,y(0)=1')例5輸出：

x=-exp(3*t)*(C1*cos(4*t)-C2*sin(4*t))

y=exp(3*t)*(C1*sin(4*t)+C2*cos(4*t))x=exp(3*t)*sin(4*t)

y=exp(3*t)*cos(4*t)MATLAB軟件求解解輸入命令：

[x,y,z]=dsolve('Dx=2*x-3*y+3*z',...'Dy=4*x-5*y+3*z',...'Dz=4*x-4*y+2*z','t');x=simple(x)%將x簡(jiǎn)化

y=simple(y)z=simple(z)結(jié)果為：x=C3*exp(2*t)+exp(-t)*C1

y=C2*exp(-2*t)+C3*exp(2*t)+exp(-t)*C1

z=C2*exp(-2*t)+C3*exp(2*t)微分方程的數(shù)值解（一）常微分方程數(shù)值解的定義

在生產(chǎn)和科研中所處理的微分方程往往很復(fù)雜且大多得不出一般解。而在實(shí)際上對(duì)初值問題，一般是要求得到解在若干個(gè)點(diǎn)上滿足規(guī)定精確度的近似值，或者得到一個(gè)滿足精確度要求的便于計(jì)算的表達(dá)式。因此，研究常微分方程的數(shù)值解法是十分必要的。返回（二）建立數(shù)值解法的一些途徑1、用差商代替導(dǎo)數(shù)

若步長(zhǎng)h較小，則有故有公式：此即歐拉法。2、使用數(shù)值積分對(duì)方程y’=f(x,y),兩邊由xi到xi+1積分，并利用梯形公式，有：實(shí)際應(yīng)用時(shí)，與歐拉公式結(jié)合使用：此即改進(jìn)的歐拉法。故有公式：3、使用泰勒公式

以此方法為基礎(chǔ)，有龍格-庫(kù)塔法、線性多步法等方法。4、數(shù)值公式的精度當(dāng)一個(gè)數(shù)值公式的截?cái)嗾`差可表示為O（hk+1）時(shí)（k為正整數(shù)，h為步長(zhǎng)），稱它是一個(gè)k階公式。k越大，則數(shù)值公式的精度越高。歐拉法是一階公式，改進(jìn)的歐拉法是二階公式。龍格-庫(kù)塔法有二階公式和四階公式。線性多步法有四階阿達(dá)姆斯外插公式和內(nèi)插公式。返回（三）用Matlab軟件求常微分方程的數(shù)值解[t，x]=solver（’f’,ts,x0,options）ode45ode23ode113ode15sode23s由待解方程寫成的m-文件名ts=[t0,tf]，t0、tf為自變量的初值和終值函數(shù)的初值ode23：組合的2/3階龍格-庫(kù)塔-芬爾格算法ode45：運(yùn)用組合的4/5階龍格-庫(kù)塔-芬爾格算法自變量值函數(shù)值用于設(shè)定誤差限(缺省時(shí)設(shè)定相對(duì)誤差10-3,絕對(duì)誤差10-6),命令為：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分別為設(shè)定的相對(duì)誤差和絕對(duì)誤差.

1、在解n個(gè)未知函數(shù)的方程組時(shí)，x0和x均為n維向量，m-文件中的待解方程組應(yīng)以x的分量形式寫成。

2、使用Matlab軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組。注意:選擇一組狀態(tài)變量注意1、建立M文件函數(shù)

functionxdot=fun(t,x,y)

xdot=[x2(t)；x3(t)；…；f(t,x1(t),x2(t),…xn(t))]；2、數(shù)值計(jì)算（執(zhí)行以下命令）

[t,x1,x2,…,xn]=ode45(‘fun',[t0,tf],[x1(0),x2(0),…,xn(0)])解:令y1=x，y2=y1’=x’1、建立m-文件vdp1000.m如下：

functiondy=vdp1000(t,y)

dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)^2)*y(2)-y(1);

2、取t0=0，tf=3000，輸入命令：

[T,Y]=ode15s('vdp1000',[03000],[20]);plot(T,Y(:,1),'-')3、結(jié)果如圖解

1、建立m-文件rigid.m如下：

functiondy=rigid(t,y)

dy=zeros(3,1);dy(1)=y(2)*y(3);dy(2)=-y(1)*y(3);dy(3)=-0.51*y(1)*y(2);2、取t0=0，tf=12，輸入命令：

[T,Y]=ode45('rigid',[012],[011]);plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')3、結(jié)果如圖圖中，y1的圖形為實(shí)線，y2的圖形為“*”線，y3的圖形為“+”線.例9Vander

pol

方程：令y1=x(t),y2=x’(t)

該方程無解析解?。?）編寫M文件(文件名為vdpol.m):

functiondy=vdpol(t,y);

dy=zeros(2,1);dy(1)=y(2);dy(2)=(1-y(1)^2)*y(2)-y(1);%或dy=[y(2);(1-y(1)^2)*y(2)-y(1)];（2）編寫程序如下:（vdj.m）

[t,y]=ode23('vdpol',[0,20],[3,0]);y1=y(:,1);%原方程的解

y2=y(:,2);plot(t,y1,t,y2,'--'

)%y1(t),y2(t)曲線圖

pause,plot(y1,y2),grid%相軌跡圖，即y2(y1)曲線

藍(lán)色曲線

——y（1）；（原方程解）

紅色曲線

——y（2）；計(jì)算結(jié)果例10考慮Lorenz模型：其中參數(shù)β=8/3，σ=10，ρ=28解：1）編寫M函數(shù)文件(lorenz.m);2)數(shù)值求解并畫三維空間的相平面軌線；

(ltest.m)1、lorenz.mfunctionxdot=lorenz(t,x)xdot=[-8/3,0,x(2);0,-10,10;-x(2),28,-1]*x;2、ltest.mx0=[000.1]';[t,x]=ode45('lorenz',[0,10],x0);plot(t,x(:,1),'-',t,x(:,2),'*',t,x(:,3),'+')pauseplot3(x(:,1),x(:,2),x(:,3)),gridon圖中，x1的圖形為實(shí)線(藍(lán)），x2的圖形為“*”線(綠)，x3的圖形為“+”線（紅）。取[t0，tf]=[0，10]。計(jì)算結(jié)果如下圖：曲線呈震蕩發(fā)散狀三維圖形的混沌狀若自變量區(qū)間取[0，20]、[0，40]，計(jì)算結(jié)果如下：觀察結(jié)果：

1、該曲線包含兩個(gè)“圓盤”，每一個(gè)都是由螺線形軌道構(gòu)成。某些軌道幾乎是垂直地離開圓盤中一個(gè)而進(jìn)入另一個(gè)。