用 Matlab 求解微分方程.ppt

《用 Matlab 求解微分方程.ppt》由會(huì)員分享，可在線閱讀，更多相關(guān)《用 Matlab 求解微分方程.ppt（28頁珍藏版）》請?jiān)谘b配圖網(wǎng)上搜索。

用 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+y^2’, ‘x’) 輸出：ans = tan(x+C1) 求特解 輸入：dsolve(‘Dy=1+y^2’, ‘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化簡 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化簡 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è)定為相對誤差 10?3，絕對誤差 10?6），程序?yàn)?options = odeset(‘reltol’, rt, ‘a(chǎn)bstol’, at) 其中rt和at分別為設(shè)定的相對誤差和絕對誤差； ⑥ 在解 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)彈頭始終對準(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)彈頭始終對準(zhǔn)乙艦，故此時(shí)直線 PQ 就是導(dǎo)彈的軌跡曲線弧 OP 在點(diǎn) P 處的切線，,,,于是有,,,,,,,,,,,,,,即,,又根據(jù)題意，弧 OP 的長度為 |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) 處擊中乙艦。,,,,,,