1、用 Matlab 求解微分方程,借助 Matlab 軟件，可以方便地求出微分方程（組）的解析解和數(shù)值解。,微分方程（組）的解析解,求微分方程（組）解析解的命令為 dsolve(eqn1, eqn2, ., x) 其中“eqni”表示第 i 個方程，“x”表示微分方程（組）中的自變量，默認時自變量為 t。此外，在“eqni”表示的方程式中，用 D 表示求微分，D2、D3 等表示求高階微分，任何 D 后所跟的字母表示因變量。,例 8.5.1 求解一階微分方程 dy/dx = 1 + y2。,求通解 輸入：dsolve(Dy=1+y2, x) 輸出：ans = tan(x+C1) 求特解 輸入：ds

2、olve(Dy=1+y2, y(0)=1, x) 輸出：ans = tan(x+1/4*pi),例 8.5.2 求解下列微分方程的通解及 y(0) = 0 和 y (0) = 15 條件下的特解,求通解 輸入：y=dsolve(D2y+4*Dy+29*y=0, x) 輸出：y = C1*exp(-2*x)*sin(5*x)+C2*exp(-2*x)*cos(5*x) 求特解 輸入：y=dsolve(D2y+4*Dy+29*y=0, y(0)=0, Dy(0)=15, x) 輸出：y = 3*exp(-2*x)*sin(5*x),例 8.5.3 求解下列微分方程組,求通解 方式一 輸入： x,y

3、,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 = C2*exp(-t)+C3*exp(2*t) y = C2*exp(-t)+C3*exp(2*t)+exp(-2*t)*C1 z = C3*exp(2*t)+exp(-2*t)*C1,方式二 輸入：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化簡 y=simple(y) z=simple(z) 輸出：x = C2/exp(t)+C3*exp(t)

4、2 y = C2*exp(-t)+C3*exp(2*t)+exp(-2*t)*C1 z = C3*exp(2*t)+exp(-2*t)*C1,求特解 輸入：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, x(0)=0, y(0)=1, z(0)=2, t)； x=simple(x) % 將x化簡 y=simple(y) z=simple(z) 輸出：x = exp(2*t)-exp(-t) y = exp(2*t)-exp(-t)+exp(-2*t) z = exp(2*t)+exp(-2*t),微分方程（組）的數(shù)值解,事

5、實上，能夠求得解析解的微分方程或微分方程組少之又少，多數(shù)情況下需要求出微分方程（組）的數(shù)值解。 Matlab中求微分方程數(shù)值解的函數(shù)有五個：ode45，ode23，ode113，ode15s，ode23s。調(diào)用格式為 t, x = solver (f, ts, x0, options),需要特別注意的是： solver 可以取以上五個函數(shù)之一，不同的函數(shù)代表不同的內(nèi)部算法：ode23 運用組合的 2/3 階龍格庫塔費爾貝算法，ode45 運用組合的 4/5 階龍格庫塔費爾貝算法。通常使用函數(shù) ode45； f 是由待解方程寫成的m文件的文件名； ts=t0, tf，t0、tf為自變量的初值和終

6、值； x0為函數(shù)的初值；, options 用于設(shè)定誤差限（可以缺省，缺省時設(shè)定為相對誤差 103，絕對誤差 106），程序為 options = odeset(reltol, rt, abstol, at) 其中rt和at分別為設(shè)定的相對誤差和絕對誤差； 在解 n 個未知函數(shù)的方程組時，x0、x 均為 n 維向量，m 文件中待解方程組應(yīng)以 x 的分量形式寫成； 使用 Matlab 軟件求數(shù)值解時，高階微分方程必須等價地變換成一階微分方程組。,例 8.5.4 求解下列微分方程,解：令 y1 = x，y2 = y1，則微分方程變?yōu)橐浑A微分方程組：,(1) 建立 m 文件 vdp1000.m 如下

7、： function dy=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,0 3000,2 0); plot(T,Y(:,1),-) 運行程序，得到如圖的結(jié)果。,例 8.5.5 求解下列微分方程組,(1) 建立 m 文件 rigid.m 如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51

8、*y(1)*y(2); (2) 取 t0=0，tf=12，輸入命令： T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+),運行程序，得到如圖的結(jié)果。圖中，y1 的圖形為實線，y2 的圖形為“*”線，y3 的圖形為“+”線。,例 8.5.6 導彈追蹤問題 設(shè)位于坐標原點的甲艦向位于 x 軸上點 A(1, 0) 處的乙艦發(fā)射導彈，導彈頭始終對準乙艦。如果乙艦以最大的速度 v0（是常數(shù)）沿平行于 y 軸的直線行駛，導彈的速度是 5v0，求導彈運行的曲線方程。又乙艦行駛多遠時，導彈將它擊中？,解：如圖所示，假設(shè)導彈在

9、t 時刻的位置為P(x(t), y(t)，乙艦位于 Q(1, v0t)。由于導彈頭始終對準乙艦，故此時直線 PQ 就是導彈的軌跡曲線弧 OP 在點 P 處的切線，,于是有,即,又根據(jù)題意，弧 OP 的長度為 |AQ| 的 5 倍，于是,消去 t，得到導彈追蹤模型如下：,下面求解這個初值問題。,解法一 解析解 利用微分方程初值問題的解析解法，得導彈的運行軌跡為：,參見下圖。,根據(jù)題意，乙艦始終沿平行于 y 軸的直線 x = 1 行駛，且由上式知，當 x = 1 時 y = 5/24，故當乙艦航行到點 (1, 5/24) 處時被導彈擊中。同時可求得被擊中時間為：t = y/v0 = 5/24v0；若 v0 = 1，則在 t = 0.21 處被擊中。,解法二 數(shù)值解 令 y1 = y，y2 = y1，將先前給出的導彈追蹤模型化為一階微分方程組,(1) 建立 m 文件 eq1.m function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(1)2)/(1-x); (2) 取 x0=0，xf=0.9999，建立主