第二講 常微分方程數(shù)值求解 —— MATLAB 求解 M

1、第二講常微分方程數(shù)值求解 MATLAB 求解2Matlab 解初值問題函數(shù) 用 Maltab自帶函數(shù) 解初值問題 求解析解：dsolve 求數(shù)值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb3符號求解dsolve符號求解4dsolve 求解析解 求解析解：dsolve y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v)其中 y 為輸出的解， eq1、eq2、. 為方程，cond1、cond2、. 為初值條件，v 為自變量。例 1：求微分方程 的通解，并驗(yàn)證。sol=dsolve(Dy+2*x*y=x*exp(

2、-x2),x)syms x diff(sol)+2*x*sol - x*exp(-x2) % 驗(yàn)證5dsolve 的使用 幾點(diǎn)說明 如果省略初值條件，則表示求通解； 如果省略自變量，則默認(rèn)自變量為 t dsolve(Dy=2*x,x); % dy/dx = 2xdsolve(Dy=2*x); % dy/dt = 2x 若找不到解析解，則提出警告，并返回空解。 微分方程中用 D 表示對 自變量 的導(dǎo)數(shù)，如：Dy y； D2y y； D3y y6dsolve 的使用 使用符號方程 導(dǎo)數(shù)：diff，如 diff(y)，diff(y,2) 等號：= 必須聲明應(yīng)變量與自變量!例 1：求微分方程 的通解，

3、并驗(yàn)證。syms y(x)sol=dsolve(diff(y)+2*x*y=x*exp(-x2)diff(y)+2*x*y - x*exp(-x2)7dsolve 舉例例 2：求微分方程 在初值條件 下的特解，并畫出解函數(shù)的圖形。sol=dsolve(x*Dy+y-exp(x)=0, y(1)=2*exp(1), x);ezplot(sol);syms y(x)sol=dsolve(diff(y)*x+y-exp(x)=0, y(1)=2*exp(1);ezplot(sol);syms y(x)sol=dsolve(diff(y)*x+y-exp(x)=0, y(1)=2*exp(sym(1)

4、;ezplot(sol);8dsolve 舉例例3：求微分方程組 在初值條件 下的特解，并畫出解函數(shù)的圖形。x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0, . x(0)=1, y(0)=0, t)ezplot(x,y,0,1.3);注：解微分方程組時(shí)，如果所給的輸出個(gè)數(shù)與方程個(gè)數(shù)相同，則方程組的解按詞典順序輸出；如果只給一個(gè)輸出，則輸出的是一個(gè)包含解的結(jié)構(gòu)（structure）類型的數(shù)據(jù)。9dsolve 舉例例：x,y=dsolve(Dx+5*x=0, Dy-3*y=0, . x(0)=1, y(0)=1, t)sol = dsolve(Dx+5*x=0, Dy-

5、3*y=0, . x(0)=1, y(0)=1, t)這里返回的 sol 是一個(gè) 結(jié)構(gòu)類型 的數(shù)據(jù)sol.x % 查看解函數(shù) x(t)sol.y % 查看解函數(shù) y(t)dsolve的輸出個(gè)數(shù)只能為一個(gè)或與方程個(gè)數(shù)相等10dsolve 舉例 使用符號方程syms x(t) y(t) sol=dsolve(diff(x)+5*x=0, diff(y)-3*y=0, . x(0)=1, y(0)=1)11dsolve 舉例dsolve(D2y=-a2*y,y(0)=1,Dy(pi/a)=0)例 4：求微分方程 的特解，初值條件為其中 a 是符號常量syms y(t) ady = diff(y);s

6、ol=dsolve(diff(y,2)=-a2*y, y(0)=1, dy(pi/a)=0)12數(shù)值求解ode45、ode23、ode113、ode23t、ode15s、ode23s、ode23tb數(shù)值求解13數(shù)值求解T, Y = solver(odefun, tspan, y0)其中 y0 為初值條件，tspan為求解區(qū)間；Matlab在數(shù)值求解時(shí)自動(dòng)對求解區(qū)間進(jìn)行分割，T (列向量) 中返回的是分割點(diǎn)的值(自變量)，Y (數(shù)組) 中返回的是這些分割點(diǎn)上的近似解，其列數(shù)等于應(yīng)變量的個(gè)數(shù)。solver 為Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、o

7、de23s、ode23t、ode23tb）沒有一種算法可以有效地解決所有的 ODE 問題，因此MATLAB 提供了多種ODE求解器，對于不同的ODE，可以調(diào)用不同的求解器。14Matlab的ODE求解器求解器ODE類型特點(diǎn)說明ode45非剛性單步法；4，5 階 R-K 方法；累計(jì)截?cái)嗾`差為 (x)3大部分場合的首選方法ode23非剛性單步法；2，3 階 R-K 方法；累計(jì)截?cái)嗾`差為 (x)3使用于精度較低的情形ode113非剛性多步法；Adams算法；高低精度均可到 10-310-6計(jì)算時(shí)間比 ode45 短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法；Gears 反向數(shù)

8、值微分；精度中等若 ode45 失效時(shí)，可嘗試使用ode23s剛性單步法；2 階Rosebrock 算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比 ode15s 短ode23tb剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短15參數(shù)說明odefun 為函數(shù)句柄，代表顯式常微分方程，可以通過匿名函數(shù)定義，或在函數(shù)文件中定義，然后通過函數(shù)句柄調(diào)用。f = (x,y) -2*y+2*x2+2*x;x,y = ode45(f, 0,0.5,1);注：也可以在 tspan 中指定對求解區(qū)間的分割，如：x,y = ode45(f, 0:0.1:0.5,1); % x=0:0.1:0.5T, Y = solver(odefun,tspan,y0) 求 的數(shù)值解，求解區(qū)間為 0,0.5例 ：16數(shù)值求解舉例如果需求解的問題是高階常微分方程，則需將其化為一階常微分方程組，此時(shí)必須用函數(shù)文件來定義該常微分方程組。令 求解 Ver der Pol 初值問題例 ：17數(shù)值求解舉例 先編寫函數(shù)文件 verderpol.mfunction dx=verderpol(t,x) % 必須是兩個(gè)輸入和一個(gè)輸出global mu; % 其它參數(shù)只能通過全局變量傳遞dx=x(2);