第八章 常微分方程的初值問題

第八章常微分方程的初值問題常微分方程（ODE）：只含有未知函數(shù)的導數(shù)的方程。ODE的階數(shù)：最高階導數(shù)的階數(shù)ODE問題初值問題：邊值問題：已知初始點處的條件已知初始點和末端點處的條件本章只討論初值問題一階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。一、求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'''例：求微分方程的通解，并驗證。>>

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：如果省略自變量，則默認自變量為t

例：

dsolve('Dy=2*x')%dy/dt=2x結(jié)果為ans=2*x*t+C1例：求微分方程在初值條件下的特解，并畫出解函數(shù)的圖形。>>

y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')結(jié)果為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')結(jié)果為ans=1/3*exp(-t)+2/3*exp(1/2*t)*cos(1/2*3^(1/2)*t)例：注4：解微分方程組時，如果所給的輸出個數(shù)與方程個數(shù)相同，則方程組的解按詞典順序輸出；如果只給一個輸出，則輸出的是一個包含解的結(jié)構(gòu)（structure）類型的數(shù)據(jù)。例：求微分方程組在初值條件下的特解[x,y]=dsolve('Dx+5*x+y=exp(t)',...'Dy-x-3*y=0','x(0)=1','y(0)=0','t')[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是一個結(jié)構(gòu)類型的數(shù)據(jù)r.x%查看解函數(shù)x(t)r.y%查看解函數(shù)y(t)dsolve的輸出個數(shù)只能為一個或與方程個數(shù)相等。例：用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ù)值解。Euler（歐拉）方法是求解一階ODE的一種簡便方法。盡管精度不高，但由于簡單，特別適用于快速編程求解。它有三種格式：二、用Euler方法求數(shù)值解（a）向前Euler法（b）改進的Euler法（c）向后Euler法介紹這些方法是為了了解初值問題求解的基本思想。一階ODE對兩邊從x0到x積分得：(積分方程）1、向前Euler法推導1：設節(jié)點為向前Euler法：用向前差分公式代替導數(shù)：此處，y(xn)表示xn處的理論解，yn表示y(xn)的近似解一階ODE對兩邊從xn

到xn+1積分得：推導2：yn近似代替y(xn)用矩形代替右邊的積分向前Euler法例求解其中C=0.27,M=70,g=9.8并畫圖C=0.27;M=70;g=9.8;x=0;y=0;h=0.1;n=0;%n為迭代次數(shù)xr(1)=t;yr(1)=v;whilex<20n=n+1;y=y+h*(-C/M*y^2+g);x=x+h;yr(n+1)=v;xr(n+1)=x;endplot(xr,yr);axis([0,20,0,60])盡管向前Euler法非常簡單，但有兩種誤差效應：第一種是截斷誤差，第二種誤差源于計算的不穩(wěn)定性。當方程的時間系數(shù)為負值，h又不是很小時，第二種誤差就會增加。帶有負時間系數(shù)的常見方程為

其精確解為該問題的向前Euler格式為：如果，數(shù)值解為正的趨于零；如果，隨n的增加，解的符號交替改變；如果，則在每一步符號變化后，解的數(shù)量級增加，最終導致結(jié)果發(fā)散。對兩邊從xn

到xn+1積分得：2、改進的Euler法yn近似代替y(xn)用梯形代替右邊的積分梯形法從n=0開始計算，每步都要求解一個關于yn+1的方程（一般是一個非線性方程），可用如下的迭代法計算：利用此算法，可得：利用（為允許誤差）來控制是否繼續(xù)迭代迭代法太麻煩，實際上，當h取得很小時，只讓上式中的第二式迭代一次就可以，即改進的Euler法（也叫歐拉預估—校正法）預估算式校正算式改進的Euler法=向前歐拉法+梯形法例2初值問題：求解的值。向前Euler法：yf(1)=10;t(1)=0;h=0.1;n=1;fprintf('tyf\n');fprintf('%3.1f%6.4f\n',t(1),yf(1));whilet(n)<1yf(n+1)=yf(n)+h*(-yf(n)^1.5+1);t(n+1)=t(n)+h;n=n+1;fprintf('%3.1f%6.4f\n',t(n),yf(n));end迭代5步結(jié)果：

tyf0.010.00000.16.93770.25.21040.34.12100.43.38440.52.8618此題無解析解ym(1)=10;t(1)=0;h=0.1;n=1;fprintf('tym\n');fprintf('%3.1f%6.4f\n',t(1),ym(1));whilet(n)<1y0(n+1)=ym(n)+h*(-ym(n)^1.5+1);ym(n+1)=ym(n)+h/2*(-y0(n+1)^1.5-ym(n)^1.5+2);t(n+1)=t(n)+h;n=n+1;fprintf('%3.1f%6.4f\n',t(n),ym(n));end改進Euler法：迭代5步結(jié)果：兩種方法的比較圖：plot(t,yf,'--',t,ym,'-');axis([0,1.2,0,10])

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

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

（b）如果解為正，則可保證數(shù)值計算結(jié)果也為正。4、二階ODE問題二階ODE的一般形式為：其中是常數(shù)或是的函數(shù)，后兩個方程為初始條件。如果與u無關，則該方程稱為線性ODE;如果三者之中有一個是u或的函數(shù)，稱為非線性的。下面著重研究用向前Euler法求解二階ODE，及MATLAB程序。所以二階ODE分解為兩個一階ODE：設：定義則上述兩個一階ODE寫為：其向前Euler法的格式為：例求解二階ODE解：設令則原方程的向量形式為：向前Euler法的格式為：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法求解5、高階ODE問題設方程:如果是常數(shù)或t的函數(shù)，則為線性的；如果至少有一個是函數(shù)，則為非線性的；初始條件：上述微分方程可轉(zhuǎn)化為四個一階ODE：向量形式為：向前Euler法的格式為：數(shù)值方法同樣適用于積分方程。如設對v微分：則原方程變?yōu)椋合蛄啃问綖椋喝?/p>