數(shù)學(xué)實(shí)驗(yàn)：求微分方程的解.ppt

1、1,實(shí)驗(yàn)四,求微分方程的解,2,自牛頓發(fā)明微積分以來(lái)，微分方程在描述事物運(yùn)動(dòng)規(guī)律上已發(fā)揮了重要的作用。實(shí)際應(yīng)用問(wèn)題通過(guò)數(shù)學(xué)建模所得到的方程，絕大多數(shù)是微分方程,由于實(shí)際應(yīng)用的需要，人們必須求解微分方程。然而能夠求得解析解的微分方程十分有限，絕大多數(shù)微分方程需要利用數(shù)值方法來(lái)近似求解,本實(shí)驗(yàn)主要研究如何用 Matlab 來(lái)計(jì)算微分方程（組）的數(shù)值解，并重點(diǎn)介紹一個(gè)求解微分方程的基本數(shù)值解法Euler折線法,問(wèn)題背景和實(shí)驗(yàn)?zāi)康?3,考慮一維經(jīng)典初值問(wèn)題,基本思想：用差商代替微商,根據(jù) Talyor 公式，y(x) 在點(diǎn) xk 處有,Euler 折線法,4,初值問(wèn)題的Euler折線法,具體步驟,等距

2、剖分,步長(zhǎng),分割求解區(qū)間,分割求解區(qū)間，差商代替微商，解代數(shù)方程,為分割點(diǎn),5,Euler 折線法舉例,例：用 Euler 法解初值問(wèn)題,取步長(zhǎng) h = (2 - 0)/n = 2/n，得差分方程,當(dāng) h=0.4，即 n=5 時(shí)，Matlab 源程序見(jiàn) fuluA.m,解,6,Euler 折線法源程序,clear; f = inline(y+2*x/y2,x,y); a = 0; b = 2; n = 5; h = (b-a)/n; x = a : h : b; y(1) = 1; for i = 1 : n y(i+1) = y(i) + h*f(x(i),y(i); end plot(x,

3、y,rs,7,Euler折線法舉例（續(xù),解析解,解析解,近似解,y=1/3*(-18-54*x+45*exp(3*x)(1/3,8,Runge-Kutta 方法,為了減小誤差，可采用以下方法,讓步長(zhǎng) h 取得更小一些,改用具有較高精度的數(shù)值方法，如,龍格-庫(kù)塔方法,Runge-Kutta (龍格-庫(kù)塔) 方法,是一類(lèi)求解常微分方程的數(shù)值方法,有多種不同的迭代格式,9,Runge-Kutta 方法,用得較多的是 四階R-K方法（教材第 98 頁(yè),其中,10,四階 R-K 方法源程序,clear f = inline(y+2*x/y2,x,y); a = 0; b = 2; n = 5; h =

4、(b-a)/n; x = a : h : b; y(1) = 1; for i = 1 : n L1 = f(x(i), y(i); L2 = f(x(i)+h/2, y(i)+L1*h/2); L3 = f(x(i)+h/2, y(i)+L2*h/2); L4 = f(x(i)+h, y(i)+L3*h); y(i+1) = y(i) + h*(L1+2*L2+2*L3+L4)/6; end fprintf(y 在 x=2 點(diǎn)的近似值為 %fn,y(n+1); plot(x,y,rs,11,Runge-Kutta 方法,12,Euler 法與 R-K法誤差比較,13,Matlab 解初值問(wèn)題

5、函數(shù),用 Maltab自帶函數(shù) 解初值問(wèn)題,求解析解：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,dsol

6、ve 的使用,幾點(diǎn)說(shuō)明,如果省略初值條件，則表示求通解,如果省略自變量，則默認(rèn)自變量為 t,dsolve(Dy=2*x,x); % dy/dx = 2x dsolve(Dy=2*x); % dy/dt = 2x,若找不到解析解，則返回其積分形式,微分方程中用 D 表示對(duì) 自變量 的導(dǎo)數(shù)，如,Dy y； D2y y； D3y y,16,dsolve 舉例,例 2：求微分方程 在初值條件 下的特解，并畫(huà)出解函數(shù)的圖形,y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y,17,dsolve 舉例,例3：求微分方程組 在初值條件 下的特解，并畫(huà)出解函數(shù)

7、的圖形,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）類(lèi)型的數(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)類(lèi)型 的數(shù)據(jù),r.x %

8、查看解函數(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í)自動(dòng)對(duì)求解區(qū)間進(jìn)行分割，T (列向量) 中返回的是分割點(diǎn)的值(自變量)，Y (數(shù)組) 中返回的是這些分割點(diǎn)上的近似解，其列數(shù)等于因變量的個(gè)數(shù)。 solver 為Matlab的ODE求解器（可以是 ode45、ode23、ode113、ode15s、ode23s、

9、ode23t、ode23tb,沒(méi)有一種算法可以有效地解決所有的 ODE 問(wèn)題，因此MATLAB 提供了多種ODE求解器，對(duì)于不同的ODE，可以調(diào)用不同的求解器,20,Matlab的ODE求解器,21,參數(shù)說(shuō)明,odefun 為顯式常微分方程，可以用命令 inline 定義，或在函數(shù)文件中定義，然后通過(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ū)間的分割，如,x,y=ode23(fun,0:0.1:0.5,1); % x=0:0.1:0.5,T,Y = solver(odefun,

10、tspan,y0,22,數(shù)值求解舉例,如果需求解的問(wèn)題是高階常微分方程，則需將其化為一階常微分方程組，此時(shí)必須用函數(shù)文件來(lái)定義該常微分方程組,令,23,數(shù)值求解舉例,先編寫(xiě)函數(shù)文件 verderpol.m,function xprime=verderpol(t,x) global mu; xprime=x(2); mu*(1-x(1)2)*x(2) - x(1,再編寫(xiě)腳本文件 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,教材 P97：練習(xí) 1 7,上機(jī)作業(yè),要求寫(xiě)實(shí)驗(yàn)報(bào)告,將所編寫(xiě)的程序分別命名為 hw41.m, hw42.m, hw43.m, hw44.m