1、求微分方程的解 數(shù)學(xué)實(shí)驗(yàn) q 自牛頓發(fā)明微積分以來(lái)，微分方程在描述事物運(yùn)自牛頓發(fā)明微積分以來(lái)，微分方程在描述事物運(yùn)動(dòng)規(guī)律上已發(fā)揮了重要的作用。實(shí)際應(yīng)用問(wèn)題通過(guò)動(dòng)規(guī)律上已發(fā)揮了重要的作用。實(shí)際應(yīng)用問(wèn)題通過(guò)數(shù)學(xué)建模所得到的方程，絕大多數(shù)是微分方程。數(shù)學(xué)建模所得到的方程，絕大多數(shù)是微分方程。q 由于實(shí)際應(yīng)用的需要，人們必須求解微分方程。由于實(shí)際應(yīng)用的需要，人們必須求解微分方程。然而能夠求得解析解的微分方程十分有限，絕大多然而能夠求得解析解的微分方程十分有限，絕大多數(shù)微分方程需要利用數(shù)值方法來(lái)近似求解。數(shù)微分方程需要利用數(shù)值方法來(lái)近似求解。q 本實(shí)驗(yàn)主要研究如何用本實(shí)驗(yàn)主要研究如何用 Matlab 來(lái)

2、計(jì)算微分方程來(lái)計(jì)算微分方程（組）的數(shù)值解，并重點(diǎn)介紹一個(gè)求解微分方程的（組）的數(shù)值解，并重點(diǎn)介紹一個(gè)求解微分方程的基本數(shù)值解法基本數(shù)值解法Euler折線法折線法。問(wèn)題背景和實(shí)驗(yàn)?zāi)康膯?wèn)題背景和實(shí)驗(yàn)?zāi)康膓 考慮一維經(jīng)典考慮一維經(jīng)典初值問(wèn)題初值問(wèn)題00 , , ( ,)() , dyf x yy xyxa bdx u 基本思想：基本思想：用差商代替微商用差商代替微商根據(jù)根據(jù) Talyor 公式，公式，y(x) 在點(diǎn)在點(diǎn) xk 處有處有2()()() ()()kkky xy xxxyxOx 1kkhxx 11()()()()( )kkkkky xy xy xy xdyO hdxhhx 21()()()

3、()kkky xy xhyxO h Euler 折線法折線法初值問(wèn)題的初值問(wèn)題的Euler折線法折線法q 具體步驟：具體步驟：等距剖分：等距剖分：0121nnaxxxxxb 步長(zhǎng)：步長(zhǎng)：1 0 1 21() /, , ,kkhxxknnba u 分割求解區(qū)間分割求解區(qū)間u 差商代替微商差商代替微商1()()kkky xy xdydxhx 1 ()()()kkky xy xh yx 得方程組：得方程組：0011 ()(,)kkkkkkyy xyyh f xyxxh 分割求解區(qū)間，差商代替微商，解代數(shù)方程分割求解區(qū)間，差商代替微商，解代數(shù)方程 為分割點(diǎn)為分割點(diǎn) 0nkkx k = 0, 1, 2,

4、 ., n-1yk 是是 y (xk) 的近似的近似Euler 折線法舉例折線法舉例例：例：用用 Euler 法解初值問(wèn)題法解初值問(wèn)題22 0 201 , ( )dyxydxyxy 取步長(zhǎng)取步長(zhǎng) h = (2 - 0)/n = 2/n，得差分方程得差分方程002110 1 2 ,(,)(/)kkkkkkkkkkxyyyh f xyyh yxyxxh 當(dāng)當(dāng) h=0.4，即即 n=5 時(shí)，時(shí)，Matlab 源程序見源程序見 Euler.m解：解：Euler 折線法源程序折線法源程序clearf=sym(y+2*x/y2);a=0; b=2;h=0.4;n=(b-a)/h; x=0; y=1;szj

5、=x,y;for i=1:n y=y+h*subs(f,x,y,x,y); x=x+h; szj=szj;x,y;endszjplot(szj(:,1),szj(:,2),or-)hold onezplot(5*exp(3*x)/3-2*x-2/3)(1/3),0,2)legend(近似解近似解,解析解析解解) Euler折線法舉例（續(xù)）折線法舉例（續(xù)）解析解：解析解：1 3352233/xyex 解析解解析解近似解近似解Runge-Kutta 方法方法q 為了減小誤差，可采用以下方法：為了減小誤差，可采用以下方法：u 讓步長(zhǎng)讓步長(zhǎng) h 取得更小一些；取得更小一些；u 改用具有較高精度的數(shù)值方

6、法：改用具有較高精度的數(shù)值方法：q 龍格龍格-庫(kù)塔方法庫(kù)塔方法Runge-Kutta (龍格龍格-庫(kù)塔庫(kù)塔) 方法方法u 是是一類一類求解常微分方程的數(shù)值方法求解常微分方程的數(shù)值方法u 有多種不同的迭代格式有多種不同的迭代格式Runge-Kutta 方法方法q 用得較多的是用得較多的是 四階四階R-K方法方法00111234 (22)/6(),kkkkyy xxxhyyh LLLL 12132432222(,)(/ ,/ )(/ ,/ )(,)kkkkkkkkLf xyLf xhyhLLf xhyhLLf xh yhL 其中其中四階四階 R-K 方法方法源程序源程序clear;f=sym(y+

7、2*x/y2);a=0; b=2; h=0.4;n=(b-a)/h;x=0; y=1; szj=x,y;for i=1:n l1=subs(f,x,y,x,y); l2=subs(f,x,y,x+h/2,y+l1*h/2); l3=subs(f,x,y,x+h/2,y+l2*h/2); l4=subs(f,x,y,x+h,y+l3*h); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=szj;x,y;endplot(szj(:,1),szj(:,2), dg-)Runge-Kutta 方法方法Euler 法與法與 R-K法誤差比較法誤差比較程序程序EuRK.mMa

8、tlab 解初值問(wèn)題解初值問(wèn)題q 用用 Maltab自帶函數(shù)自帶函數(shù) 解初值問(wèn)題解初值問(wèn)題u 求解析解：求解析解：dsolveu 求數(shù)值解：求數(shù)值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tbdsolve 求解析解求解析解q dsolve 的使用的使用y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v)其中其中 y 為輸出，為輸出， eq1、eq2、.為微分方程，為微分方程，cond1、cond2、.為初值條件，為初值條件，v 為自變量。為自變量。例例 1：求微分方程求微分方程 的通解，并驗(yàn)證。的通解，并驗(yàn)證。

9、22xdyxyxedx y=dsolve(Dy+2*x*y=x*exp(-x2),x) syms x; diff(y)+2*x*y - x*exp(-x2)dsolve 的使用的使用q 幾點(diǎn)說(shuō)明幾點(diǎn)說(shuō)明l 如果省略初值條件，則表示求通解；如果省略初值條件，則表示求通解；l 如果省略自變量，則默認(rèn)自變量為如果省略自變量，則默認(rèn)自變量為 t dsolve(Dy=2*x,x); dy/dx = 2xdsolve(Dy=2*x); dy/dt = 2xl 若找不到解析解，則返回其積分形式。若找不到解析解，則返回其積分形式。l 微分方程中用微分方程中用 D 表示對(duì)表示對(duì) 自變量自變量 的導(dǎo)數(shù)，如：的導(dǎo)數(shù)

10、，如：Dy y； D2y y； D3y ydsolve 舉例舉例例例 2：求微分方程求微分方程 在初值條件在初值條件 下的特解，并畫出解函數(shù)的圖形。下的特解，并畫出解函數(shù)的圖形。0 xxyye y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y);12( )ye dsolve 舉例舉例例例3：求微分方程組求微分方程組 在初值條件在初值條件 下的特解，并畫出解函數(shù)的圖形。下的特解，并畫出解函數(shù)的圖形。530tdxxyedtdyxydt x,y=dsolve(Dx+5*x+y=exp(t),Dy-x-3*y=0, . x(0)=1, y(0)=0,

11、 t)ezplot(x,y,0,1.3);0010|ttxy 注：解微分方程組時(shí)，如果所給的輸出個(gè)數(shù)與方程個(gè)數(shù)相同，注：解微分方程組時(shí)，如果所給的輸出個(gè)數(shù)與方程個(gè)數(shù)相同，則方程組的解則方程組的解按詞典順序按詞典順序輸出；如果只給一個(gè)輸出，則輸出輸出；如果只給一個(gè)輸出，則輸出的是一個(gè)包含解的的是一個(gè)包含解的結(jié)構(gòu)結(jié)構(gòu)（structure）類型的數(shù)據(jù)。類型的數(shù)據(jù)。v dsolve 舉例舉例例：例：x,y=dsolve(Dx+5*x=0,Dy-3*y=0, . x(0)=1, y(0)=1,t) r = dsolve(Dx+5*x=0,Dy-3*y=0, . x(0)=1, y(0)=1,t)這里返回

12、的這里返回的 r 是一個(gè)是一個(gè) 結(jié)構(gòu)類型結(jié)構(gòu)類型 的數(shù)據(jù)的數(shù)據(jù)r.x %查看解函數(shù)查看解函數(shù) x(t)r.y %查看解函數(shù)查看解函數(shù) y(t)只有很少一部分微分方程（組）能求出解析解。只有很少一部分微分方程（組）能求出解析解。大部分微分方程（組）只能利用大部分微分方程（組）只能利用數(shù)值方法數(shù)值方法求數(shù)值解。求數(shù)值解。dsolve的輸出個(gè)數(shù)只能為一個(gè)的輸出個(gè)數(shù)只能為一個(gè) 或或 與方程個(gè)數(shù)相等。與方程個(gè)數(shù)相等。Matlab函數(shù)數(shù)值求解函數(shù)數(shù)值求解T,Y = solver(odefun,tspan,y0,options)其中其中 y0 為初值條件，為初值條件，tspan為求解區(qū)間為求解區(qū)間(自變量的

13、初值和終值自變量的初值和終值)；Matlab在數(shù)值求解時(shí)在數(shù)值求解時(shí)自動(dòng)對(duì)求解區(qū)間進(jìn)行分割自動(dòng)對(duì)求解區(qū)間進(jìn)行分割，T (向量向量) 中返中返回的是分割點(diǎn)的值回的是分割點(diǎn)的值(自變量自變量)，Y (向量向量) 中返回的是解函數(shù)在這中返回的是解函數(shù)在這些分割點(diǎn)上的函數(shù)值。些分割點(diǎn)上的函數(shù)值。options設(shè)定誤差限設(shè)定誤差限(缺省缺省:相對(duì)相對(duì)10-3, 絕絕對(duì)對(duì)10-6),命令：命令：options=odeset（reltol,rt,abstol,at）, rt，at：分別為設(shè)定的相對(duì)誤差和絕對(duì)誤差分別為設(shè)定的相對(duì)誤差和絕對(duì)誤差.solver 為為Matlab的的ODE求解器求解器（可以是（可以

14、是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）沒有一種算法可以有效地解決所有的沒有一種算法可以有效地解決所有的 ODE 問(wèn)題，因此問(wèn)題，因此Matlab 提供了多種提供了多種ODE求解器求解器，對(duì)于不同的對(duì)于不同的ODE，可以調(diào)用不同的，可以調(diào)用不同的求解器求解器。ode只能直接求解一階微分，只能直接求解一階微分，高階方法：把一個(gè)高階方法：把一個(gè)N階階微分方程轉(zhuǎn)化為微分方程轉(zhuǎn)化為N個(gè)一階微分方程組。個(gè)一階微分方程組。Matlab提供的提供的ODE求解器求解器求解器求解器 ODE類型類型特點(diǎn)特點(diǎn)說(shuō)明說(shuō)明ode45非剛性非剛性通常情況下最好的單

15、步法；通常情況下最好的單步法；4，5 階階 R-K 方法；誤差為方法；誤差為 (x)3大部分場(chǎng)合的大部分場(chǎng)合的首選方首選方法法(速度慢速度慢)ode23非剛性非剛性單步法；單步法；2，3 階階 R-K 方法；方法；累計(jì)截?cái)嗾`差為累計(jì)截?cái)嗾`差為 (x)3使用于精度較低的情使用于精度較低的情形形(比比ode45更小步長(zhǎng)更小步長(zhǎng))ode113非剛性非剛性多步法；多步法；Adams算法；高低精度算法；高低精度均可到均可到 10-310-6計(jì)算時(shí)間比計(jì)算時(shí)間比 ode45 短短,適合連續(xù)的系統(tǒng)適合連續(xù)的系統(tǒng)ode23t適度剛性適度剛性 采用梯形算法采用梯形算法適度剛性情形適度剛性情形ode15s剛性剛性

16、多步法；多步法；Gears 反向數(shù)值微分；反向數(shù)值微分；精度中等（精度中等（可用可用隱式隱式）若若 ode45 失效時(shí)，失效時(shí)，可嘗試使用可嘗試使用(速度快速度快)ode23s剛性剛性單步法；單步法；2 階階Rosebrock 算法；算法；低精度低精度當(dāng)精度較低時(shí)，計(jì)算當(dāng)精度較低時(shí)，計(jì)算時(shí)間比時(shí)間比 ode15s 短短ode23tb剛性剛性梯形算法；低精度梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算當(dāng)精度較低時(shí)，計(jì)算時(shí)間比時(shí)間比ode15s短短參數(shù)說(shuō)明參數(shù)說(shuō)明odefun 為為顯式常微分方程顯式常微分方程，可以用命令，可以用命令 inline 定義，或定義，或在在函數(shù)文件函數(shù)文件中定義，然后通過(guò)函數(shù)句柄

17、調(diào)用。中定義，然后通過(guò)函數(shù)句柄調(diào)用。fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,0,0.5,1);注：注：也可以在也可以在 tspan 中指定對(duì)求解區(qū)間的分割，如：中指定對(duì)求解區(qū)間的分割，如：x,y=ode23(fun,0:0.1:0.5,1); %此時(shí)此時(shí) x=0:0.1:0.5T,Y = solver(odefun,tspan,y0) 求初值問(wèn)題求初值問(wèn)題 的數(shù)值解，求解范的數(shù)值解，求解范圍為圍為 0,0.5222201( )dyyxxdxy 例例 4：數(shù)值求解舉例數(shù)值求解舉例如果需求解的問(wèn)題是如果需求解的問(wèn)題是高階高階常微分方程，則需將其化為常

18、微分方程，則需將其化為一階常一階常微分方程組微分方程組，此時(shí)需用，此時(shí)需用函數(shù)文件函數(shù)文件來(lái)定義該常微分方程組。來(lái)定義該常微分方程組。122212112101 00 7/()( ),( ),dxdtxdxdtxxxxx 令令 ，則原方程可化為，則原方程可化為12,dyxy xdt 求解求解 Ver der Pol 初值問(wèn)題初值問(wèn)題2221001 00 7()( ),( ),d ydyyydtdtyy 例例 5：數(shù)值求解舉例數(shù)值求解舉例l 先編寫函數(shù)文件先編寫函數(shù)文件 verderpol.mfunction xprime=verderpol(t,x)%t,x分別給出的時(shí)間向量和狀態(tài)向量分別給出的

19、時(shí)間向量和狀態(tài)向量global mu;xprime=x(2); mu*(1-x(1)2)*x(2) - x(1);l 再編寫腳本文件再編寫腳本文件 vdpl.m，在命令窗口直接運(yùn)行該文件。，在命令窗口直接運(yùn)行該文件。 clear;global mu;mu=7;y0=1;0;t,x=ode45(verderpol,0,40,y0); plot(t,x(:,1),rd);hold ont1,x1=ode23(verderpol,0,40,y0); plot(t1,x1(:,1), bo);l 先編寫函數(shù)文件先編寫函數(shù)文件 verderpol1.ml 再編寫腳本文件再編寫腳本文件 vdpl1.m，在

20、命令窗口直接運(yùn)行該文件。，在命令窗口直接運(yùn)行該文件。 clear;global mu;mu=7;y0=1;0;t,x=ode45(verderpol1,0,40,y0); plot(t,x(:,1),rd);hold ont1,x1=ode15s(verderpol1,0,40,y0); plot(t1,x1(:,1), bo);function dy=verderpol1(t,y) global mu; dy=zeros(2,1); dy(1)=y(2); dy(2)=mu*(1-y(1)2)*y(2)-y(1);數(shù)值求解舉例數(shù)值求解舉例例例6 6：求解微分方程：求解微分方程, 1)0(,

21、1ytyy要求求解析解和數(shù)值解，并進(jìn)行比較。要求求解析解和數(shù)值解，并進(jìn)行比較。s=dsolve(Dy=-y+t+1,y(0)=1,t)%解析解解析解%數(shù)值解數(shù)值解，先編寫，先編寫M文件文件fun1.mfunction f=fun1(t,y)f=-y+t+1;%m文件中輸入文件中輸入clear;t=0:0.1:1;y=t+exp(-t); plot(t,y);hold on %保留已經(jīng)畫好的圖形，兩個(gè)圖形合并在一起保留已經(jīng)畫好的圖形，兩個(gè)圖形合并在一起t,y=ode45(fun1,0,1,1); plot(t,y,r*) %畫數(shù)值解圖形，用畫數(shù)值解圖形，用*畫畫xlabel(t)ylabel(y)Matlab 求解微分方程小結(jié)求解微分方程小結(jié)q Matlab 函數(shù)函數(shù)u 求解析解（求解析解（通解通解或特解），用或特解），用 dsolveu 求數(shù)值解（特解），用求數(shù)值解（特解），用 ode45、ode23 .q Matlab 編程編程u Euler 折線法折線法u Runga-Kutta 方法方法例例 洛倫茲模型洛倫茲模型由如下常微分方程組描述由如下常微分方程組描述( (蝴蝶效應(yīng)蝴蝶效應(yīng)) zyxydtdzzydtdyyzxdtdx )(取取 =8/3， =10， =28。初值初值：x(0)=0，y(0