用 Matlab 求解微分方程PPT課件

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

2、e(Dy=1+y2, x)輸出：ans = tan(x+C1)求特解輸入：dsolve(Dy=1+y2, y(0)=1, x)輸出：ans = tan(x+1/4*pi)2021/3/94例 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)*s

3、in(5*x)029422ydxdydxyd2021/3/95例 8.5.3 求解下列微分方程組 zyxdtdzzyxdtdyzyxdtdx2443543322021/3/96求通解 方式一 輸入： 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)*C12021/3/97方式二輸入：x,y,z=dsolve(Dx=2*x-3*y

4、+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)*C12021/3/98求特解 輸入：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化

5、簡 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)2021/3/99微分方程（組）的數(shù)值解微分方程（組）的數(shù)值解 事實上，能夠求得解析解的微分方程或微分方程組少之又少，多數(shù)情況下需要求出微分方程（組）的數(shù)值解。 Matlab中求微分方程數(shù)值解的函數(shù)有五個：ode45，ode23，ode113，ode15s，ode23s。調(diào)用格式為t, x = solver (f, ts, x0, options)2021/3/910 需要特別注意的是： s

6、olver 可以取以上五個函數(shù)之一，不同的函數(shù)代表不同的內(nèi)部算法：ode23 運用組合的 2/3 階龍格庫塔費爾貝算法，ode45 運用組合的 4/5 階龍格庫塔費爾貝算法。通常使用通常使用函數(shù)函數(shù) ode45； f 是由待解方程寫成的m文件的文件名； ts=t0, tf，t0、tf為自變量的初值和終值； x0為函數(shù)的初值；2021/3/911 options 用于設(shè)定誤差限（可以缺省，缺省時設(shè)定為相對誤差 103，絕對誤差 106），程序為options = odeset(reltol, rt, abstol, at)其中rt和at分別為設(shè)定的相對誤差和絕對誤差； 在解 n 個未知函數(shù)的方程

7、組時，x0、x 均為 n 維向量，m 文件中待解方程組應(yīng)以 x 的分量形式寫成； 使用 Matlab 軟件求數(shù)值解時，高階微分方程必須等價地變換成一階微分方程組。 2021/3/9122021/3/913例 8.5.4 求解下列微分方程0)0( ; 2)0(0)1 (1000222xxxdtdxxdtxd解：令 y1 = x，y2 = y1，則微分方程變?yōu)橐浑A微分方程組： 0)0(, 2)0()1 (1000211221221yyyyyyyy2021/3/914 (1) 建立 m 文件 vdp1000.m 如下： function dy=vdp1000(t,y) dy=zeros(2,1);

8、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é)果。 2021/3/9152021/3/916例 8.5.5 求解下列微分方程組1)0(, 1)0(, 0)0(51. 0321213312321yyyyyyyyyyyy2021/3/917(1) 建立 m 文件 rigid.m 如下： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3)

9、; 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),+)2021/3/918 運行程序，得到如圖的結(jié)果。圖中，y1 的圖形為實線，y2 的圖形為“*”線，y3 的圖形為“+”線。2021/3/919例 8.5.6 導(dǎo)彈追蹤問題導(dǎo)彈追蹤問題 設(shè)位于坐標原點的甲艦向位于 x 軸上點 A(1, 0) 處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對準乙艦。如果乙艦以最大的速度 v0（是常數(shù)）沿平行于 y 軸的直線

10、行駛，導(dǎo)彈的速度是 5v0，求導(dǎo)彈運行的曲線方程。又乙艦行駛多遠時，導(dǎo)彈將它擊中？2021/3/920解：如圖所示，假設(shè)導(dǎo)彈在 t 時刻的位置為P(x(t), y(t)，乙艦位于 Q(1, v0t)。由于導(dǎo)彈頭始終對準乙艦，故此時直線 PQ 就是導(dǎo)彈的軌跡曲線弧 OP 在點 P 處的切線， 2021/3/921于是有 xytvy10即yyxtv)1 (0又根據(jù)題意，弧 OP 的長度為 |AQ| 的 5 倍，于是 tvdxyx002512021/3/922消去 t，得到導(dǎo)彈追蹤模型如下： 下面求解這個初值問題。 0)0( , 0)0(151)1 (2yyyyx2021/3/923解法一 解析解

11、利用微分方程初值問題的解析解法，得導(dǎo)彈的運行軌跡為： 參見下圖。245)1 (125)1 (855654xxy2021/3/92400.10.20.30.40.50.60.70.80.9100.050.10.150.20.252021/3/925 根據(jù)題意，乙艦始終沿平行于 y 軸的直線 x = 1 行駛，且由上式知，當 x = 1 時 y = 5/24，故當乙艦航行到點 (1, 5/24) 處時被導(dǎo)彈擊中。同時可求得被擊中時間為：t = y/v0 = 5/24v0；若 v0 = 1，則在 t = 0.21 處被擊中。 2021/3/926解法二 數(shù)值解 令 y1 = y，y2 = y1，將先前給出的導(dǎo)彈追蹤模型化為一階微分方程組 0)0( , 0)0()1/(1512121221yyxyyyy2021/3/927(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),