常微分方程的初值問(wèn)題

常微分方程的初值問(wèn)題第一頁(yè)，共五十三頁(yè)，2022年，8月28日一階ODE問(wèn)題的形式1、如果f與y無(wú)關(guān)，則計(jì)算變?yōu)榈?章（數(shù)值積分）中討論的任一種直接積分方法應(yīng)用

初始條件2、如果f是y的函數(shù)，積分過(guò)程將不同于前者。若f是y

的線(xiàn)性函數(shù)，如：f=ay+b

其中a,b是常數(shù)或是

t的函數(shù)，此時(shí)原方程稱(chēng)為線(xiàn)性O(shè)DE若f不是線(xiàn)性函數(shù)，方程就稱(chēng)為非線(xiàn)性O(shè)DE。第二頁(yè)，共五十三頁(yè)，2022年，8月28日一、求ODE的解析解dsolve[輸出變量列表]=dsolve(‘eq1’,‘eq2’,...,‘eqn’,‘cond1’,‘cond2’,...,‘condn’,‘v1,v2,…vn')其中

eq1、eq2、...、eqn為微分方程，cond1、cond2、...、condn為初值條件，v1,v2,…,vn

為自變量。注1：微分方程中用

D表示對(duì)自變量的導(dǎo)數(shù)，如：Dyy'；

D2yy''；

D3yy'''第三頁(yè)，共五十三頁(yè)，2022年，8月28日例：求微分方程的通解，并驗(yàn)證。>>

y=dsolve('Dy+2*x*y=x*exp(-x^2)','x')

結(jié)果為y=(1/2*x^2+C1)*exp(-x^2)>>

symsx;diff(y)+2*x*y-x*exp(-x^2)

結(jié)果為ans=0注2：如果省略初值條件，則表示求通解；注3：如果省略自變量，則默認(rèn)自變量為t

例：

dsolve('Dy=2*x')%dy/dt=2x結(jié)果為ans=2*x*t+C1第四頁(yè)，共五十三頁(yè)，2022年，8月28日例：求微分方程在初值條件下的特解，并畫(huà)出解函數(shù)的圖形。>>

y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')結(jié)果為y=(exp(x)+exp(1))/x

>>ezplot(y)%此時(shí)的y已經(jīng)是符號(hào)變量，故不用ezplot(‘y’)dsolve('D3y=-y','y(0)=1,Dy(0)=0,D2y(0)=0')結(jié)果為ans=1/3*exp(-t)+2/3*exp(1/2*t)*cos(1/2*3^(1/2)*t)例：第五頁(yè)，共五十三頁(yè)，2022年，8月28日注4：解微分方程組時(shí)，如果所給的輸出個(gè)數(shù)與方程個(gè)數(shù)相同，則方程組的解按詞典順序輸出；如果只給一個(gè)輸出，則輸出的是一個(gè)包含解的結(jié)構(gòu)（structure）類(lèi)型的數(shù)據(jù)。例：求微分方程組在初值條件下的特解[x,y]=dsolve('Dx+5*x+y=exp(t)',...'Dy-x-3*y=0','x(0)=1','y(0)=0','t')第六頁(yè)，共五十三頁(yè)，2022年，8月28日[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0',...'x(0)=1','y(0)=0','t')或r=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0',...'x(0)=1','y(0)=0','t')這里返回的r是一個(gè)結(jié)構(gòu)類(lèi)型的數(shù)據(jù)r.x%查看解函數(shù)

x(t)r.y%查看解函數(shù)

y(t)dsolve的輸出個(gè)數(shù)只能為一個(gè)或與方程個(gè)數(shù)相等。第七頁(yè)，共五十三頁(yè)，2022年，8月28日例：用dsolve求解著名的VanderPol方程

>>symsmu;>>y=dsolve('D2y+mu*(y^2-1)*Dy+y')結(jié)果：Warning,cannotfindanexplicitsolutiony=&where(_a,{[diff(_b(_a),_a)*_b(_a)+mu*_b(_a)*_a^2-mu*_b(_a)+_a=0,_b(_a)=diff(y(t),t),_a=y(t),y(t)=_a,t=Int(1/_b(_a),_a)+C1]})注5：若找不到解析解，則返回其積分形式。只有很少一部分微分方程（組）能求出解析解。

大部分微分方程（組）只能利用數(shù)值方法求數(shù)值解。第八頁(yè)，共五十三頁(yè)，2022年，8月28日

Euler（歐拉）方法是求解一階ODE的一種簡(jiǎn)便方法。盡管精度不高，但由于簡(jiǎn)單，特別適用于快速編程求解。它有三種格式：二、用Euler方法求數(shù)值解（a）向前Euler法（b）改進(jìn)的Euler法（c）向后Euler法介紹這些方法是為了了解初值問(wèn)題求解的基本思想。第九頁(yè)，共五十三頁(yè)，2022年，8月28日一階ODE對(duì)兩邊從x0到x

積分得：

(積分方程）第十頁(yè)，共五十三頁(yè)，2022年，8月28日1、向前Euler法推導(dǎo)1：設(shè)節(jié)點(diǎn)為向前Euler法：用向前差分公式代替導(dǎo)數(shù)：此處，y(xn)表示xn處的理論解，yn表示y(xn)的近似解第十一頁(yè)，共五十三頁(yè)，2022年，8月28日一階ODE對(duì)兩邊從xn

到xn+1

積分得：推導(dǎo)2：yn近似代替y(xn)用矩形代替右邊的積分向前Euler法第十二頁(yè)，共五十三頁(yè)，2022年，8月28日例求出解析解和數(shù)值解并畫(huà)圖比較先用dsolve求解析解y=dsolve('Dy=y+2*x/(y^2)','y(0)=1','x')結(jié)果為y=1/3*(-18-54*x+45*exp(3*x))^(1/3)解析解：第十三頁(yè)，共五十三頁(yè)，2022年，8月28日f(shuō)unction[x,y]=Euler_bf(fun,x0,y0,xmax,h)%fun為顯式一階微分方程右端的函數(shù)%x0,y0為初始條件,滿(mǎn)足y(x0)=y0%xmax為x的取值最大值%h為將x等分后的步長(zhǎng)N=(xmax-x0)/h+1;%N為總的節(jié)點(diǎn)數(shù)x(1)=x0;y(1)=y0;forn=1:N-1x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n));end下面再用向前歐拉法數(shù)值求解，為此先編寫(xiě)向前歐拉法的程序第十四頁(yè)，共五十三頁(yè)，2022年，8月28日最后通過(guò)圖形比較用向前歐拉得到的數(shù)值解和解析解的誤差

fun=inline('y+2*x/y^2','x','y');[x1,y1]=Euler_bf(fun,0,1,2,0.1);[x2,y2]=Euler_bf(fun,0,1,2,0.05);x=0:0.01:2;y=1/3*(-18-54*x+45*exp(3*x)).^(1/3);plot(x1,y1,'*',x2,y2,'x',x,y)axis([0,2,0,9])第十五頁(yè)，共五十三頁(yè)，2022年，8月28日向前歐拉方法截?cái)嗾`差為第十六頁(yè)，共五十三頁(yè)，2022年，8月28日對(duì)兩邊從xn

到xn+1

積分得：2、改進(jìn)的Euler法yn近似代替y(xn)用梯形代替右邊的積分梯形法第十七頁(yè)，共五十三頁(yè)，2022年，8月28日從n=0開(kāi)始計(jì)算，每步都要求解一個(gè)關(guān)于yn+1的方程（一般是一個(gè)非線(xiàn)性方程），可用如下的迭代法計(jì)算：利用此算法，可得：利用（為允許誤差）來(lái)控制是否繼續(xù)迭代第十八頁(yè)，共五十三頁(yè)，2022年，8月28日迭代法太麻煩，實(shí)際上，當(dāng)h取得很小時(shí)，只讓上式中的第二式迭代一次就可以，即改進(jìn)的Euler法（也叫歐拉預(yù)估—校正法）預(yù)估算式校正算式改進(jìn)的Euler法=向前歐拉法+梯形法第十九頁(yè)，共五十三頁(yè)，2022年，8月28日向后Euler法依賴(lài)于向后差分近似，其格式為：

3、向后Euler法精度與向前歐拉法相同。如果f為非線(xiàn)性函數(shù)，則與改進(jìn)的Euler算法一樣，在每一步中需要采用迭代法。該方法有兩個(gè)優(yōu)點(diǎn)：

（a）絕對(duì)穩(wěn)定；

（b）如果解為正，則可保證數(shù)值計(jì)算結(jié)果也為正。第二十頁(yè)，共五十三頁(yè)，2022年，8月28日三、龍格－庫(kù)塔（Runge-kutta）方法

Euler法的一個(gè)重要缺陷是精度太低。為了獲得高精度，就要減小h，這不僅會(huì)增加計(jì)算時(shí)間，也會(huì)產(chǎn)生舍入誤差。1、二階Runge-kutta方法第二十一頁(yè)，共五十三頁(yè)，2022年，8月28日或其實(shí)就是歐拉預(yù)估—校正方法（或改進(jìn)的歐拉法）第二十二頁(yè)，共五十三頁(yè)，2022年，8月28日例用改進(jìn)的歐拉法（即二階龍格-庫(kù)塔法）數(shù)值求解并與向前歐拉法、解析解畫(huà)圖比較前面已求出解析解：第二十三頁(yè)，共五十三頁(yè)，2022年，8月28日f(shuō)unction[x,y]=Euler_mend(fun,x0,y0,xmax,h)%fun為顯式一階微分方程右端的函數(shù)%x0,y0為初始條件,滿(mǎn)足y(x0)=y0%xmax為x的取值最大值%h為將x等分后的步長(zhǎng)N=(xmax-x0)/h+1;%N為總的節(jié)點(diǎn)數(shù)x(1)=x0;y(1)=y0;forn=1:N-1x(n+1)=x(n)+h;k1=h*feval(fun,x(n),y(n));k2=h*feval(fun,x(n+1),y(n)+k1);y(n+1)=y(n)+1/2*(k1+k2);end先編寫(xiě)改進(jìn)的歐拉法的程序：第二十四頁(yè)，共五十三頁(yè)，2022年，8月28日再分別調(diào)用Euler_bf.m和Euler_mend.m求解：

fun=inline('y+2*x/y^2','x','y');[x1,y1]=Euler_bf(fun,0,1,2,0.1);[x2,y2]=Euler_mend(fun,0,1,2,0.1);x=0:0.01:2;y=1/3*(-18-54*x+45*exp(3*x)).^(1/3);plot(x1,y1,'*',x2,y2,‘s',x,y)axis([0,2,0,9])第二十五頁(yè)，共五十三頁(yè)，2022年，8月28日第二十六頁(yè)，共五十三頁(yè)，2022年，8月28日3、三階龍格－庫(kù)塔方法常見(jiàn)的三階龍格－庫(kù)塔方法的格式為：二階龍格－庫(kù)塔方法截?cái)嗾`差為精度不高，希望通過(guò)增加計(jì)算f的次數(shù)提高截?cái)嗾`差的階數(shù)，為此引入第二十七頁(yè)，共五十三頁(yè)，2022年，8月28日4、四階龍格－庫(kù)塔方法常見(jiàn)的四階龍格-庫(kù)塔方法有兩種，一種基于1/3辛普森法，格式：第二十八頁(yè)，共五十三頁(yè)，2022年，8月28日另一種基于3/8辛普森法，格式：第二十九頁(yè)，共五十三頁(yè)，2022年，8月28日例用二階龍格-庫(kù)塔法和四階龍格-庫(kù)塔法數(shù)值求解并與、解析解畫(huà)圖比較前面已求出解析解：先來(lái)編寫(xiě)四階龍格-庫(kù)塔法（基于1/3辛普森法）的程序：第三十頁(yè)，共五十三頁(yè)，2022年，8月28日f(shuō)unction[x,y]=RK4(fun,x0,y0,xmax,h)%fun為顯式一階微分方程右端的函數(shù)%x0,y0為初始條件,滿(mǎn)足y(x0)=y0%xmax為x的取值最大值%h為將x等分后的步長(zhǎng)N=(xmax-x0)/h+1;%N為總的節(jié)點(diǎn)數(shù)x(1)=x0;y(1)=y0;forn=1:N-1x(n+1)=x(n)+h;k1=h*feval(fun,x(n),y(n));k2=h*feval(fun,x(n)+1/2*h,y(n)+k1/2);k3=h*feval(fun,x(n)+1/2*h,y(n)+k2/2);k4=h*feval(fun,x(n+1),y(n)+k3);y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4);end第三十一頁(yè)，共五十三頁(yè)，2022年，8月28日f(shuō)un=inline('y+2*x/y^2','x','y');[x1,y1]=Euler_mend(fun,0,1,2,0.1);[x2,y2]=RK4(fun,0,1,2,0.1);x=0:0.1:2;y=1/3*(-18-54*x+45*exp(3*x)).^(1/3);plot(x1,y1,'*',x2,y2,'s',x,y)axis([0,2,0,9])再分別調(diào)用Euler_mend.m和RK4.m求解：第三十二頁(yè)，共五十三頁(yè)，2022年，8月28日從圖形上看，似乎二階龍格—庫(kù)塔法與四階龍格-庫(kù)塔法同樣接近解析解，故再?gòu)臄?shù)值結(jié)果比較看看哪種誤差小第三十三頁(yè)，共五十三頁(yè)，2022年，8月28日err=[abs(y'-y1'),abs(y'-y2')]結(jié)果為err=000.00090.00000.00160.00000.00220.00000.00260.00000.00310.00000.00350.00000.00400.00000.00470.00000.00540.00000.00620.00000.00730.00000.00840.00000.00980.00000.01150.00000.01330.00000.01550.00000.01810.00000.02100.00000.02430.00000.02810.0000第三十四頁(yè)，共五十三頁(yè)，2022年，8月28日四、二階ODE問(wèn)題二階ODE的一般形式為：其中是常數(shù)或是的函數(shù)，后兩個(gè)方程為初始條件。如果與u無(wú)關(guān)，則該方程稱(chēng)為線(xiàn)性O(shè)DE;如果三者之中有一個(gè)是u或的函數(shù)，稱(chēng)為非線(xiàn)性的。下面著重研究用向前Euler法求解二階ODE，及MATLAB程序。第三十五頁(yè)，共五十三頁(yè)，2022年，8月28日所以二階ODE分解為兩個(gè)一階ODE：設(shè)：第三十六頁(yè)，共五十三頁(yè)，2022年，8月28日定義則上述兩個(gè)一階ODE寫(xiě)為：其向前Euler法的格式為：第三十七頁(yè)，共五十三頁(yè)，2022年，8月28日例求解二階ODE解：設(shè)令則原方程的向量形式為：向前Euler法的格式為：第三十八頁(yè)，共五十三頁(yè)，2022年，8月28日h=0.05;t_max=5;n=1;y(:,1)=[0;1];t(1)=0;whilet(n)<t_maxy(:,n+1)=y(:,n)+h*f_def(y(:,n),t);t(n+1)=t(n)+h;n=n+1;endaxis([05-11]);plot(t,y(1,:),t,y(2,:),':')functionf=f_def(y,t)f=[y(2);(-5*abs(y(2))*y(2)-20*y(1))];先用函數(shù)文件定義f(u,v,t)再用向前Euler法求解第三十九頁(yè)，共五十三頁(yè)，2022年，8月28日第四十頁(yè)，共五十三頁(yè)，2022年，8月28日五、matlab相關(guān)命令數(shù)值求解ODE[X,Y]

=求解函數(shù)(fun,[x0,xf],y0,option，p1,p2,….)其中：（1）fun是用inline或函數(shù)文件定義的顯式常微分方程的函數(shù)名函數(shù)文件格式：或inline格式：functionyd=函數(shù)名(x,y,flag,p1,p2,…)yd=…….函數(shù)名=inline(‘顯式方程’,’x’,’y’,’flag’,’p1’,’p2’,….)x為自變量，y為因變量，yd為y的導(dǎo)數(shù),flag用于控制求解過(guò)程，p1，p2為方程中的參數(shù)第四十一頁(yè)，共五十三頁(yè)，2022年，8月28日[X,Y]

=求解函數(shù)(fun,[x0,xf],y0,option，p1,p2,….)其中：（2）[x0,xf]為求解區(qū)間，y0為初值條件；（3）option（可省略）為控制選項(xiàng)（用于設(shè)置誤差限，求解方程最大允許步長(zhǎng)等等），（4）p1,p2為微分方程中的附加參數(shù)（5）Matlab在數(shù)值求解時(shí)自動(dòng)對(duì)求解區(qū)間進(jìn)行分割，X(向量)中返回的是分割點(diǎn)的值(自變量)，Y(向量)中返回的是解函數(shù)在這些分割點(diǎn)上的函數(shù)值。（6）求解函數(shù)可以是ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb第四十二頁(yè)，共五十三頁(yè)，2022年，8月28日求解函數(shù)ODE類(lèi)型特點(diǎn) 說(shuō)明ode45非剛性單步法；4，5階R-K方法；累計(jì)截?cái)嗾`差為(△x)3大部分場(chǎng)合的首選方法ode23非剛性單步法；2，3階R-K方法；累計(jì)截?cái)嗾`差為(△x)3使用于精度較低的情形ode113非剛性多步法；Adams算法；高低精度均可到10-3～10-6

計(jì)算時(shí)間比ode45

短ode23t適度剛性采用梯形算法適度剛性情形ode15s剛性多步法；Gear’s反向數(shù)值微分；精度中等若ode45

失效時(shí)，可嘗試使用ode23s剛性單步法；2階Rosebrock算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短ode23tb剛性梯形算法；低精度當(dāng)精度較低時(shí)，計(jì)算時(shí)間比ode15s短沒(méi)有一種算法可以有效地解決所有的ODE問(wèn)題，因此MATLAB提供了多種ODE求解函數(shù)，對(duì)于不同的ODE，可以調(diào)用不同的求解函數(shù)。第四十三頁(yè)，共五十三頁(yè)，2022年，8月28日求初值問(wèn)題的數(shù)值解，求解范圍為[0,0.5]例先定義需要求解的顯式微分方程為一個(gè)函數(shù)functionyd=fun1(x,y)yd=-2*y+2*x^2+2*x最后在命令窗口調(diào)用函數(shù)求解[x,y]=ode23(‘fun1’,[0,0.5],1);第四十四頁(yè)，共五十三頁(yè)，2022年，8月28日f(shuō)un2=inline('-2*y+2*x^2+2*x','x','y');求初值問(wèn)題的數(shù)值解，求解范圍為[0,0.5]例先定義需要求解的顯式微分方程為一個(gè)函數(shù)在命令窗口用inline定義最后在命令窗口調(diào)用函數(shù)求解[x,y]=ode23(fun2,[0,0.5],1);第四十五頁(yè)，共五十三頁(yè)，2022年，8月28日如果需求解的問(wèn)題是高階常微分方程，則需將其化為一階常微分方程組，此時(shí)需用函數(shù)文件來(lái)定義該常微分方程組。令

，則原方程可化為求解VerderPol初值問(wèn)題例：前面演示過(guò)，該方程沒(méi)有解析解，該方程不是一階顯式微分方程，故需要進(jìn)行變換第四十六頁(yè)，共五十三頁(yè)，2022年，8月28日先用函數(shù)文件定義一階顯式微分方程組functiony=vdp_eq1(t,x,flag,mu)y=[x(2);-mu*(x(1)^2-1)*x(2)-x(1)];再編寫(xiě)腳本文件

vdpl.m，在命令窗口直接運(yùn)行該文件。x0=[-0.2;-0.7];tf=20;mu=1;[t1,y1]=ode45('vdp_eq1',[0,tf],x0,[],mu);mu=2;[t2,y2]=ode45('vdp_eq1',[0,tf],x0,[],mu);plot(t1,y1,t2,y2,':')figure;plot(y1(:,1),y1(:,2),y2(:,1),y2(:,2),':')法1：第四十七頁(yè)，共五十三頁(yè)，2022年，8月28日vdp_eq2=inline('[x(2);-mu*(x(1)^2-1)*x(2)-x(1)]',…'t','x','flag','mu');x0=[-0.2;-0.7];tf=20;mu=1;[t1,y1]=ode45(vdp_eq2,[0,tf],x0,[],mu);mu=2