常微分方程的初值問題

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

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

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

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

t的函數(shù)，此時原方程稱為線性ODE若f不是線性函數(shù)，方程就稱為非線性ODE。第二頁，共五十三頁，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表示對自變量的導數(shù)，如：Dyy'；

D2yy''；

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

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

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

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

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

例：

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

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

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

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

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

>>symsmu;>>y=dsolve('D2y+mu*(y^2-1)*Dy+y')結果：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ù)值解。第八頁，共五十三頁，2022年，8月28日

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

積分得：

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

到xn+1

積分得：推導2：yn近似代替y(xn)用矩形代替右邊的積分向前Euler法第十二頁，共五十三頁，2022年，8月28日例求出解析解和數(shù)值解并畫圖比較先用dsolve求解析解y=dsolve('Dy=y+2*x/(y^2)','y(0)=1','x')結果為y=1/3*(-18-54*x+45*exp(3*x))^(1/3)解析解：第十三頁，共五十三頁，2022年，8月28日function[x,y]=Euler_bf(fun,x0,y0,xmax,h)%fun為顯式一階微分方程右端的函數(shù)%x0,y0為初始條件,滿足y(x0)=y0%xmax為x的取值最大值%h為將x等分后的步長N=(xmax-x0)/h+1;%N為總的節(jié)點數(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ù)值求解，為此先編寫向前歐拉法的程序第十四頁，共五十三頁，2022年，8月28日最后通過圖形比較用向前歐拉得到的數(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])第十五頁，共五十三頁，2022年，8月28日向前歐拉方法截斷誤差為第十六頁，共五十三頁，2022年，8月28日對兩邊從xn

到xn+1

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

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

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

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

Euler法的一個重要缺陷是精度太低。為了獲得高精度，就要減小h，這不僅會增加計算時間，也會產(chǎn)生舍入誤差。1、二階Runge-kutta方法第二十一頁，共五十三頁，2022年，8月28日或其實就是歐拉預估—校正方法（或改進的歐拉法）第二十二頁，共五十三頁，2022年，8月28日例用改進的歐拉法（即二階龍格-庫塔法）數(shù)值求解并與向前歐拉法、解析解畫圖比較前面已求出解析解：第二十三頁，共五十三頁，2022年，8月28日function[x,y]=Euler_mend(fun,x0,y0,xmax,h)%fun為顯式一階微分方程右端的函數(shù)%x0,y0為初始條件,滿足y(x0)=y0%xmax為x的取值最大值%h為將x等分后的步長N=(xmax-x0)/h+1;%N為總的節(jié)點數(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先編寫改進的歐拉法的程序：第二十四頁，共五十三頁，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])第二十五頁，共五十三頁，2022年，8月28日第二十六頁，共五十三頁，2022年，8月28日3、三階龍格－庫塔方法常見的三階龍格－庫塔方法的格式為：二階龍格－庫塔方法截斷誤差為精度不高，希望通過增加計算f的次數(shù)提高截斷誤差的階數(shù)，為此引入第二十七頁，共五十三頁，2022年，8月28日4、四階龍格－庫塔方法常見的四階龍格-庫塔方法有兩種，一種基于1/3辛普森法，格式：第二十八頁，共五十三頁，2022年，8月28日另一種基于3/8辛普森法，格式：第二十九頁，共五十三頁，2022年，8月28日例用二階龍格-庫塔法和四階龍格-庫塔法數(shù)值求解并與、解析解畫圖比較前面已求出解析解：先來編寫四階龍格-庫塔法（基于1/3辛普森法）的程序：第三十頁，共五十三頁，2022年，8月28日function[x,y]=RK4(fun,x0,y0,xmax,h)%fun為顯式一階微分方程右端的函數(shù)%x0,y0為初始條件,滿足y(x0)=y0%xmax為x的取值最大值%h為將x等分后的步長N=(xmax-x0)/h+1;%N為總的節(jié)點數(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第三十一頁，共五十三頁，2022年，8月28日fun=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求解：第三十二頁，共五十三頁，2022年，8月28日從圖形上看，似乎二階龍格—庫塔法與四階龍格-庫塔法同樣接近解析解，故再從數(shù)值結果比較看看哪種誤差小第三十三頁，共五十三頁，2022年，8月28日err=[abs(y'-y1'),abs(y'-y2')]結果為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第三十四頁，共五十三頁，2022年，8月28日四、二階ODE問題二階ODE的一般形式為：其中是常數(shù)或是的函數(shù)，后兩個方程為初始條件。如果與u無關，則該方程稱為線性ODE;如果三者之中有一個是u或的函數(shù)，稱為非線性的。下面著重研究用向前Euler法求解二階ODE，及MATLAB程序。第三十五頁，共五十三頁，2022年，8月28日所以二階ODE分解為兩個一階ODE：設：第三十六頁，共五十三頁，2022年，8月28日定義則上述兩個一階ODE寫為：其向前Euler法的格式為：第三十七頁，共五十三頁，2022年，8月28日例求解二階ODE解：設令則原方程的向量形式為：向前Euler法的格式為：第三十八頁，共五十三頁，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法求解第三十九頁，共五十三頁，2022年，8月28日第四十頁，共五十三頁，2022年，8月28日五、matlab相關命令數(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的導數(shù),flag用于控制求解過程，p1，p2為方程中的參數(shù)第四十一頁，共五十三頁，2022年，8月28日[X,Y]

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

計算時間比ode45

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

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

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

vdpl.m，在命令窗口直接運行該文件。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：第四十七頁，共五十三頁，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