求微分方程的解

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

2、atlab 來計算微分方程來計算微分方程（組）的數(shù)值解，并重點介紹一個求解微分方程的（組）的數(shù)值解，并重點介紹一個求解微分方程的基本數(shù)值解法基本數(shù)值解法euler折線法折線法。問題背景和實驗目的問題背景和實驗目的3q 考慮一維經(jīng)典考慮一維經(jīng)典初值問題初值問題00 , , ( ,)() , dyf x yy xyxa bdx u 基本思想：基本思想：用差商代替微商用差商代替微商根據(jù)根據(jù) talyor 公式，公式，y(x) 在點在點 xk 處有處有2()()() ()()kkky xy xxxyxox 1kkhxx 11()()()()( )kkkkky xy xy xy xdyo hdxhhx

3、21()()()()kkky xy xhyxo h euler 折線法折線法4初值問題的初值問題的euler折線法折線法q 具體步驟：具體步驟：等距剖分：等距剖分：0121nnaxxxxxb 步長：步長：1() /, 0,1,2,1kkkhxxbann u 分割求解區(qū)間分割求解區(qū)間分割求解區(qū)間，差商代替微商，解代數(shù)方程分割求解區(qū)間，差商代替微商，解代數(shù)方程 為分割點為分割點 0nkkx u 差商代替微商差商代替微商1()()kkky xy xdydxhx 1 ()()()kkky xy xh yx 得迭代格式：得迭代格式：0011() (,) kkkkkkyy xyyh f xyxxh k =

4、 0, 1, 2, ., n-1yk 是是 y (xk) 的近似的近似5euler 折線法舉例折線法舉例例：例：用用 euler 法解初值問題法解初值問題22 0 201 , ( )dyxydxyxy 取步長取步長 h = (2 - 0)/n = 2/n，得差分方程得差分方程00110 1 2 ,(,)(/)kkkkkkkkkkxyyyh f xyyh yxyxxh 當當 h=0.4，即即 n=5 時，時，matlab 源程序見源程序見 fulua.m解：解：6euler 折線法源程序折線法源程序clearf=sym(y+2*x/y2);a=0; b=2;h=0.4;n=(b-a)/h+1;

5、% n=(b-a)/h;x=0; y=1;szj=x,y;for i=1:n-1 % 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-)7 euler折線法舉例（續(xù)）折線法舉例（續(xù)）解析解：解析解：解析解解析解近似解近似解y=1/3*(-18-54*x+45*exp(3*x)(1/3)8runge-kutta 方法方法q 為了減小誤差，可采用以下方法：為了減小誤差，可采用以下方法：u 讓步長讓步長 h 取得更小一些；取得更小一些；u 改用具有較高精度的數(shù)值方法：改用具有較高精度的數(shù)值

6、方法：q 龍格龍格-庫塔方法庫塔方法runge-kutta (龍格龍格-庫塔庫塔) 方法方法u 是是一類一類求解常微分方程的數(shù)值方法求解常微分方程的數(shù)值方法u 有多種不同的迭代格式有多種不同的迭代格式9runge-kutta 方法方法q 用得較多的是用得較多的是 四階四階r-k方法方法00111234 (22)/6(),kkkkyy xxxhyyh llll 12132432222(,)(/ ,/ )(/ ,/ )(,)kkkkkkkklf xylf xhyhllf xhyhllf xh yhl 其中其中10四階四階 r-k 方法源程序方法源程序clear;f=sym(y+2*x/y2);a=

7、0; b=2; h=0.4;n=(b-a)/h+1; % n=(b-a)/h;x=0; y=1; szj=x,y;for i=1:n-1 % 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-)11runge-kutta 方法方法12euler 法與法與 r-k法誤差

8、比較法誤差比較13matlab 解初值問題解初值問題q 用用 maltab自帶函數(shù)自帶函數(shù) 解初值問題解初值問題u 求解析解：求解析解：dsolveu 求數(shù)值解：求數(shù)值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb14dsolve 求解析解求解析解q dsolve 的使用的使用y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v)其中其中 y 為輸出，為輸出， eq1、eq2、.為微分方程，為微分方程，cond1、cond2、.為初值條件，為初值條件，v 為自變量。為自變量。例例 1：求微分方程求微分方程 的通解

9、，并驗證。的通解，并驗證。22xdyxyxedx y=dsolve(dy+2*x*y=x*exp(-x2),x) syms x; diff(y)+2*x*y - x*exp(-x2)15dsolve 的使用的使用q 幾點說明幾點說明l 如果省略初值條件，則表示求通解；如果省略初值條件，則表示求通解；l 如果省略自變量，則默認自變量為如果省略自變量，則默認自變量為 t dsolve(dy=2*x,x); dy/dx = 2xdsolve(dy=2*x); dy/dt = 2xl 若找不到解析解，則返回其積分形式。若找不到解析解，則返回其積分形式。l 微分方程中用微分方程中用 d 表示對表示對 自

10、變量自變量 的導數(shù)，如：的導數(shù)，如：dy y； d2y y； d3y y16dsolve 舉例舉例例例 2：求微分方程求微分方程 在初值條件在初值條件 下的特解，并畫出解函數(shù)的圖形。下的特解，并畫出解函數(shù)的圖形。0 xxyye y=dsolve(x*dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y);12( )ye 17dsolve 舉例舉例例例3：求微分方程組求微分方程組 在初值條件在初值條件 下的特解，并畫出解函數(shù)的圖形。下的特解，并畫出解函數(shù)的圖形。530tdxxyedtdyxydt x,y=dsolve(dx+5*x+y=exp(t),dy-x-3*y=0

11、, . x(0)=1, y(0)=0, t)ezplot(x,y,0,1.3);0010|ttxy 注：解微分方程組時，如果所給的輸出個數(shù)與方程個數(shù)相同，注：解微分方程組時，如果所給的輸出個數(shù)與方程個數(shù)相同，則方程組的解則方程組的解按詞典順序按詞典順序輸出；如果只給一個輸出，則輸出輸出；如果只給一個輸出，則輸出的是一個包含解的的是一個包含解的結構結構（structure）類型的數(shù)據(jù)。類型的數(shù)據(jù)。18 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

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

13、atlab在數(shù)值求解在數(shù)值求解時時自動對求解區(qū)間進行分割自動對求解區(qū)間進行分割，t (列向量列向量) 中返回的是分割點中返回的是分割點的值的值(自變量自變量)，y (數(shù)組數(shù)組) 中返回的是這些分割點上的近似解，中返回的是這些分割點上的近似解，其列數(shù)等于因變量的個數(shù)。其列數(shù)等于因變量的個數(shù)。solver 為為matlab的的ode求解器求解器（可以是（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）沒有一種算法可以有效地解決所有的沒有一種算法可以有效地解決所有的 ode 問題，因此問題，因此matlab 提供了多種提供了多種ode求解器求解

14、器，對于不同的對于不同的ode，可以調(diào)用不同的可以調(diào)用不同的求解器求解器。20matlab的的ode求解器求解器求解器求解器 ode類型類型特點特點說明說明ode45非剛性非剛性單步法；單步法；4 4，5 5 階階 r-k r-k 方法；方法；累計截斷誤差為累計截斷誤差為 (x)3大部分場合的大部分場合的首選方法首選方法ode23非剛性非剛性單步法；單步法；2 2，3 3 階階 r-k r-k 方法；方法；累計截斷誤差為累計截斷誤差為 (x)3使用于精度較低的情形使用于精度較低的情形ode113非剛性非剛性多步法；多步法；adams算法；高低精算法；高低精度均可到度均可到 10-310-6計算

15、時間比計算時間比 ode45 短短ode23t適度剛性適度剛性 采用梯形算法采用梯形算法適度剛性情形適度剛性情形ode15s剛性剛性多步法；多步法；gears gears 反向數(shù)值微反向數(shù)值微分；精度中等分；精度中等若若 ode45 失效時，可失效時，可嘗試使用嘗試使用ode23s剛性剛性單步法；單步法；2 2 階階rosebrockrosebrock 算法；算法；低精度低精度當精度較低時，計算時當精度較低時，計算時間比間比 ode15s 短短ode23tb剛性剛性梯形算法；低精度梯形算法；低精度當精度較低時，計算時當精度較低時，計算時間比間比ode15s短短21參數(shù)說明參數(shù)說明odefun

16、為為顯式常微分方程顯式常微分方程，可以用命令，可以用命令 inline 定義，或定義，或在在函數(shù)文件函數(shù)文件中定義，然后通過函數(shù)句柄調(diào)用。中定義，然后通過函數(shù)句柄調(diào)用。fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,0,0.5,1);注：注：也可以在也可以在 tspan 中指定對求解區(qū)間的分割，如：中指定對求解區(qū)間的分割，如：x,y=ode23(fun,0:0.1:0.5,1); %此時此時 x=0:0.1:0.5t,y = solver(odefun,tspan,y0) 求求 的數(shù)值解，求解區(qū)間為的數(shù)值解，求解區(qū)間為 0,0.5222201( )dyy

17、xxdxy 例例 ：22數(shù)值求解舉例數(shù)值求解舉例如果需求解的問題是如果需求解的問題是高階高階常微分方程，則需將其化為常微分方程，則需將其化為一階常一階常微分方程組微分方程組，此時必須用，此時必須用函數(shù)文件函數(shù)文件來定義該常微分方程組。來定義該常微分方程組。122212112(1)7(0)1, (0)0, dxxdtdxxxxdtxx 令令 12 xydyxdt 求解求解 ver der pol 初值問題初值問題222(1)0 7, (0)1, (0)0d ydyyydtdtyy 例例 ：23數(shù)值求解舉例數(shù)值求解舉例l 先編寫函數(shù)文件先編寫函數(shù)文件 verderpol.mfunction xpr

18、ime=verderpol(t,x)global mu;xprime=x(2); mu*(1-x(1)2)*x(2) - x(1);l 再編寫腳本文件再編寫腳本文件 vdpl.m，在命令窗口直接運行該文件，在命令窗口直接運行該文件 clear;global mu;mu=7; y0=1; 0;%t,x=ode45(verderpol, 0,40, y0); t,x=ode45(verderpol, 0,40, y0);plot(t,x(:,1);24求解微分方程小結求解微分方程小結q matlab 函數(shù)函數(shù)u 求解析解（求解析解（通解通解或特解），用或特解），用 dsolveu 求數(shù)值解（特解）

19、，用求數(shù)值解（特解），用 ode45、ode23 .q matlab 編程編程u euler 折線法折線法u runga-kutta 方法方法25如何定義函數(shù)文件如何定義函數(shù)文件t,y = ode45(odefun,tspan,y0)t,y = ode23(odefun,tspan,y0)當當 是函數(shù)是函數(shù)向量時呢？向量時呢？微分方程組微分方程組odefun那么那么odefun就是就是26如何定義函數(shù)文件如何定義函數(shù)文件t,y = ode45(odefun,tspan,y0)t,y = ode23(odefun,tspan,y0)ode45、ode23 等函數(shù)可用于求解等函數(shù)可用于求解顯式常微

20、分方程顯式常微分方程當當 是向量函數(shù)時，所對應的方程即為是向量函數(shù)時，所對應的方程即為微分方程組微分方程組odefun27舉例說明舉例說明fun=inline(-2*y+2*x2+2*x,x,y);x,y=ode23(fun,0,0.5,1);例：求初值問題例：求初值問題 的數(shù)值解。的數(shù)值解。解法一：解法一：使用使用 inline 定義微分方程定義微分方程 odefunodefun 為方程右端項為方程右端項 f(t,y)可以用可以用 inline 定義（只適合于定義（只適合于單個方程單個方程的情形的情形）通過通過函數(shù)文件函數(shù)文件定義，然后用函數(shù)句柄調(diào)用（適合所有情形）定義，然后用函數(shù)句柄調(diào)用（

21、適合所有情形）注：自變量必須在前面，因變量在后面！注：自變量必須在前面，因變量在后面！28舉例說明（單個方程）舉例說明（單個方程）function dy = myfun1(x,y) dy = -2*y+2*x2+2*x;解法二：解法二：通過通過函數(shù)文件函數(shù)文件定義微分方程定義微分方程 odefun1、先編寫函數(shù)文件、先編寫函數(shù)文件 myfun1.mclear;x,y=ode23(myfun1,0,0.5,1);2、編寫主文件編寫主文件 main1.m或直接在或直接在 matlab 命令窗口輸入上面的語句。命令窗口輸入上面的語句。29舉例說明（方程組）舉例說明（方程組）解解：此時只能通過：此時只

22、能通過函數(shù)文件函數(shù)文件定義微分方程定義微分方程 odefun例：求例：求 , , , , 的數(shù)值解。的數(shù)值解。function dy = myfun2(t,y) dy = zeros(3,1); % dy must be a column vector! dy(1) = y(2) * y(3); dy(2) = -y(1) * y(3); dy(3) = -0.51 * y(1) * y(2);1、先編寫函數(shù)文件、先編寫函數(shù)文件 myfun2.mclear;t,y=ode45(myfun2,0,12,0,1,1);2、編寫主文件編寫主文件 main2.mdy = y(2)*y(3); -y(1

23、)*y(3); . -0.51*y(1)*y(2);30思考思考function dy = myfun2(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);1、函數(shù)文件、函數(shù)文件 myfun2.m能不能寫成下面形式？能不能寫成下面形式？function dy = myfun2(t,x,y,z) dy = zeros(3,1); dy(1) = y * z; dy(2) = -x * z; dy(3) = -0.51 * x * y;x31說明說明odefun變量屬性必須變量屬性必須一一對應！一一對應！function dy = myfun2(t,y)如果是常微分方程組，如果是常微分方程組，y 就是列向量！就是列向量！dy 必須是必須是列列向量向量，長度為方程組的個數(shù)，通常與，長度為方程組的個數(shù)，通常與y的長度相同！的長度相同！函數(shù)中的輸入?yún)?shù)和輸出參數(shù)是函數(shù)中的輸入?yún)?shù)和輸出參數(shù)是形參形參，