1,實(shí)驗(yàn)四 求微分方程的解,數(shù)學(xué)實(shí)驗(yàn),2,自牛頓發(fā)明微積分以來，微分方程在描述事物運(yùn)動規(guī)律上已發(fā)揮了重要的作用。實(shí)際應(yīng)用問題通過數(shù)學(xué)建模所得到的方程，絕大多數(shù)是微分方程。,由于實(shí)際應(yīng)用的需要，人們必須求解微分方程。然而能夠求得解析解的微分方程十分有限，絕大多數(shù)微分方程需要利用數(shù)值方法來近似求解。,本實(shí)驗(yàn)主要研究如何用 Matlab 來計(jì)算微分方程（組）的數(shù)值解，并重點(diǎn)介紹一個(gè)求解微分方程的基本數(shù)值解法Euler折線法。,問題背景和實(shí)驗(yàn)?zāi)康?3,考慮一維經(jīng)典初值問題,基本思想：用差商代替微商,根據(jù) Talyor 公式，y(x) 在點(diǎn) xk 處有,Euler 折線法,4,初值問題的Euler折線法,具體步驟：,等距剖分：,步長：,分割求解區(qū)間,分割求解區(qū)間，差商代替微商，解代數(shù)方程,為分割點(diǎn),5,Euler 折線法舉例,例：用 Euler 法解初值問題,取步長 h = (2 - 0)/n = 2/n，得差分方程,當(dāng) h=0.4，即 n=5 時(shí)，Matlab 源程序見 fuluA.m,解：,6,Euler 折線法源程序,clear f=sym(y+2*x/y2); a=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 y=y+h*subs(f,x,y,x,y); x=x+h; szj=szj;x,y; end szj plot(szj(:,1),szj(:,2),or-),7,Euler折線法舉例（續(xù)）,解析解：,解析解,近似解,y=1/3*(-18-54*x+45*exp(3*x)(1/3),8,Runge-Kutta 方法,為了減小誤差，可采用以下方法：,讓步長 h 取得更小一些；,改用具有較高精度的數(shù)值方法：,龍格-庫塔方法,Runge-Kutta (龍格-庫塔) 方法,是一類求解常微分方程的數(shù)值方法,有多種不同的迭代格式,9,Runge-Kutta 方法,用得較多的是 四階R-K方法,其中,10,四階 R-K 方法源程序,clear; f=sym(y+2*x/y2); a=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; end plot(szj(:,1),szj(:,2), dg-),11,Runge-Kutta 方法,12,Euler 法與 R-K法誤差比較,13,Matlab 解初值問題,用 Maltab自帶函數(shù) 解初值問題,求解析解：dsolve,求數(shù)值解： ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb,14,dsolve 求解析解,dsolve 的使用,y=dsolve(eq1,eq2, . ,cond1,cond2, . ,v),其中 y 為輸出， eq1、eq2、.為微分方程，cond1、cond2、.為初值條件，v 為自變量。,例 1：求微分方程 的通解，并驗(yàn)證。, y=dsolve(Dy+2*x*y=x*exp(-x2),x), syms x; diff(y)+2*x*y - x*exp(-x2),15,dsolve 的使用,幾點(diǎn)說明,如果省略初值條件，則表示求通解；,如果省略自變量，則默認(rèn)自變量為 t,dsolve(Dy=2*x,x); dy/dx = 2x dsolve(Dy=2*x); dy/dt = 2x,若找不到解析解，則返回其積分形式。,微分方程中用 D 表示對 自變量 的導(dǎo)數(shù)，如：,Dy y； D2y y； D3y y,16,dsolve 舉例,例 2：求微分方程 在初值條件 下的特解，并畫出解函數(shù)的圖形。, y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y);,17,dsolve 舉例,例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ù)。,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(0)=1, y(0)=1,t),這里返回的 r 是一個(gè) 結(jié)構(gòu)類型 的數(shù)據(jù),r.x %查看解函數(shù) x(t) r.y %查看解函數(shù) y(t),只有很少一部分微分方程（組）能求出解析解。 大部分微分方程（組）只能利用數(shù)值方法求數(shù)值解。,dsolve的輸出個(gè)數(shù)只能為一個(gè) 或 與方程個(gè)數(shù)相等,19,數(shù)值求解,T,Y = solver(odefun,tspan,y0),其中 y0 為初值條件，tspan為求解區(qū)間；Matlab在數(shù)值求解時(shí)自動對求解區(qū)間進(jìn)行分割，T (列向量) 中返回的是分割點(diǎn)的值(自變量)，Y (數(shù)組) 中返回的是這些分割點(diǎn)上的近似解，其列數(shù)等于因變量的個(gè)數(shù)。 solver 為Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）,沒有一種算法可以有效地解決所有的 ODE 問題，因此MATLAB 提供了多種ODE求解器，對于不同的ODE，可以調(diào)用不同的求解器。,20,Matlab的ODE求解器,21,參數(shù)說明,odefun 為顯式常微分方程，可以用命令 inline 定義，或在函數(shù)文件中定義，然后通過函數(shù)句柄調(diào)用。,fun=inline(-2*y+2*x2+2*x,x,y); x,y=ode23(fun,0,0.5,1);,注：也可以在 tspan 中指定對求解區(qū)間的分割，如：,x,y=ode23(fun,0:0.1:0.5,1); %此時(shí) x=0:0.1:0.5,T,Y = solver(odefun,tspan,y0),22,數(shù)值求解舉例,如果需求解的問題是高階常微分方程，則需將其化為一階常微分方程組，此時(shí)必須用函數(shù)文件來定義該常微分方程組。,令,23,數(shù)值求解舉例,先編寫函數(shù)文件 verderpol.m,function xprime=verderpol(t,x) global mu; xprime=x(2); mu*(1-x(1)2)*x(2) - x(1);,再編寫腳本文件 vdpl.m，在命令窗口直接運(yùn)行該文件,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,求解微分方程小結(jié),Matlab 函數(shù),求解析解（通解或特解），用 dsolve 求數(shù)值解（特解），用 ode45、ode23 .,Matlab 編程,Euler 折線法 Runga-Kutta 方法,25,如何定義函數(shù)文件,T,Y = ode45(odefun,tspan,y0) T,Y = ode23(odefun,tspan,y0),當(dāng) 是函數(shù)向量時(shí)呢？,那么odefun就是,26,如何定義函數(shù)文件,T,Y = ode45(odefun,tspan,y0) T,Y = ode23(odefun,tspan,y0),ode45、ode23 等函數(shù)可用于求解顯式常微分方程 當(dāng) 是向量函數(shù)時(shí)，所對應(yīng)的方程即為微分方程組,odefun,27,舉例說明,fun=inline(-2*y+2*x2+2*x,x,y); x,y=ode23(fun,0,0.5,1);,例：求初值問題 的數(shù)值解。,解法一：使用 inline 定義微分方程 odefun,odefun 為方程右端項(xiàng) f(t,y) 可以用 inline 定義（只適合于單個(gè)方程的情形） 通過函數(shù)文件定義，然后用函數(shù)句柄調(diào)用（適合所有情形）,注：自變量必須在前面，因變量在后面！,28,舉例說明（單個(gè)方程）,function dy = myfun1(x,y) dy = -2*y+2*x2+2*x;,解法二：通過函數(shù)文件定義微分方程 odefun,1、先編寫函數(shù)文件 myfun1.m,clear; x,y=ode23(myfun1,0,0.5,1);,2、編寫主文件 main1.m,或直接在 Matlab 命令窗口輸入上面的語句。,29,舉例說明（方程組）,解：此時(shí)只能通過函數(shù)文件定義微分方程 odefun,例：求 , , 的數(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ù)文件 myfun2.m,clear; T,Y=ode45(myfun2,0,12,0,1,1);,2、編寫主文件 main2.m,dy = y(2)*y(3); -y(1)*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ù)文件 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;,X,31,說明,odefun,變量屬性必須一一對應(yīng)！,function dy = myfun2(t,y),如果是常微分方程組，y 就是列向量！,dy 必須是列向量，長度為方程組的個(gè)數(shù)，通常與y的長度相同！,函數(shù)中的輸入?yún)?shù)和輸出參數(shù)是形參，名字可以任意取，但必須滿足上述條件。即輸入?yún)?shù)有兩個(gè)，第一個(gè)表示自變量，第二個(gè)是由因變量組成的列向量，輸出參數(shù)必須是列向量。,32,例,function dy = myfun3(t,y) dy = zeros(2,1); dy(1) = y(2) + t; dy(2) = t - 2;,例：解初值問題： ， ，,function out = myfun3(t,y) out = y(2) + t; t 2;,function yprime = myfun3(x,y) yprime = y(2) + x; x 2;,clear; T,Y=ode45(myfun3,0,10,1,1);,2、主文件 main3.m,1、函數(shù)文件 myfun3.m,33,高階常微分方程,高階常微分方程,例： Van der Pol 初值問題,令 ，則原方程可化為,一階常微分方程組,變量替換 化為,參數(shù)怎么處理？,用