用 Matlab 求解微分方程.ppt

用 Matlab 求解微分方程,借助 Matlab 軟件，可以方便地求出微分方程（組）的解析解和數(shù)值解。,微分方程（組）的解析解,求微分方程（組）解析解的命令為 dsolve(eqn1, eqn2, ., x) 其中“eqni”表示第 i 個(gè)方程，“x”表示微分方程（組）中的自變量，默認(rèn)時(shí)自變量為 t。此外，在“eqni”表示的方程式中，用 D 表示求微分，D2、D3 等表示求高階微分，任何 D 后所跟的字母表示因變量。,例 8.5.1 求解一階微分方程 dy/dx = 1 + y2。,求通解 輸入：dsolve(Dy=1+y2, x) 輸出：ans = tan(x+C1) 求特解 輸入：dsolve(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,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化簡(jiǎn) y=simple(y) z=simple(z) 輸出：x = C2/exp(t)+C3*exp(t)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化簡(jiǎn) 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ù)值解,事實(shí)上，能夠求得解析解的微分方程或微分方程組少之又少，多數(shù)情況下需要求出微分方程（組）的數(shù)值解。 Matlab中求微分方程數(shù)值解的函數(shù)有五個(gè)：ode45，ode23，ode113，ode15s，ode23s。調(diào)用格式為 t, x = solver (f, ts, x0, options),需要特別注意的是： solver 可以取以上五個(gè)函數(shù)之一，不同的函數(shù)代表不同的內(nèi)部算法：ode23 運(yùn)用組合的 2/3 階龍格庫塔費(fèi)爾貝算法，ode45 運(yùn)用組合的 4/5 階龍格庫塔費(fèi)爾貝算法。通常使用函數(shù) ode45； f 是由待解方程寫成的m文件的文件名； ts=t0, tf，t0、tf為自變量的初值和終值； x0為函數(shù)的初值；, options 用于設(shè)定誤差限（可以缺省，缺省時(shí)設(shè)定為相對(duì)誤差 103，絕對(duì)誤差 106），程序?yàn)?options = odeset(reltol, rt, abstol, at) 其中rt和at分別為設(shè)定的相對(duì)誤差和絕對(duì)誤差； 在解 n 個(gè)未知函數(shù)的方程組時(shí)，x0、x 均為 n 維向量，m 文件中待解方程組應(yīng)以 x 的分量形式寫成； 使用 Matlab 軟件求數(shù)值解時(shí)，高階微分方程必須等價(jià)地變換成一階微分方程組。,例 8.5.4 求解下列微分方程,解：令 y1 = x，y2 = y1，則微分方程變?yōu)橐浑A微分方程組：,(1) 建立 m 文件 vdp1000.m 如下： 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),'-') 運(yùn)行程序，得到如圖的結(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*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),'+'),運(yùn)行程序，得到如圖的結(jié)果。圖中，y1 的圖形為實(shí)線，y2 的圖形為“*”線，y3 的圖形為“+”線。,例 8.5.6 導(dǎo)彈追蹤問題 設(shè)位于坐標(biāo)原點(diǎn)的甲艦向位于 x 軸上點(diǎn) A(1, 0) 處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對(duì)準(zhǔn)乙艦。如果乙艦以最大的速度 v0（是常數(shù)）沿平行于 y 軸的直線行駛，導(dǎo)彈的速度是 5v0，求導(dǎo)彈運(yùn)行的曲線方程。又乙艦行駛多遠(yuǎn)時(shí)，導(dǎo)彈將它擊中？,解：如圖所示，假設(shè)導(dǎo)彈在 t 時(shí)刻的位置為P(x(t), y(t)，乙艦位于 Q(1, v0t)。由于導(dǎo)彈頭始終對(duì)準(zhǔn)乙艦，故此時(shí)直線 PQ 就是導(dǎo)彈的軌跡曲線弧 OP 在點(diǎn) P 處的切線，,于是有,即,又根據(jù)題意，弧 OP 的長(zhǎng)度為 |AQ| 的 5 倍，于是,消去 t，得到導(dǎo)彈追蹤模型如下：,下面求解這個(gè)初值問題。,解法一 解析解 利用微分方程初值問題的解析解法，得導(dǎo)彈的運(yùn)行軌跡為：,參見下圖。,根據(jù)題意，乙艦始終沿平行于 y 軸的直線 x = 1 行駛，且由上式知，當(dāng) x = 1 時(shí) y = 5/24，故當(dāng)乙艦航行到點(diǎn) (1, 5/24) 處時(shí)被導(dǎo)彈擊中。同時(shí)可求得被擊中時(shí)間為：t = y/v0 = 5/24v0；若 v0 = 1，則在 t = 0.21 處被擊中。,解法二 數(shù)值解 令 y1 = y，y2 = y1，將先前給出的導(dǎo)彈追蹤模型化為一階微分方程組,(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，建立主程序 ff6.m 如下： x0=0, xf=0.9999 x,y=ode45('eq1',x0 xf,0 0); plot(x,y(:,1), 'b.') hold on y=0:0.01:2; plot(1,y, 'b*'),運(yùn)行程序，得到如圖所示的結(jié)果。從而得出結(jié)論：導(dǎo)彈大致在 (1, 0.2) 處擊中乙艦。,