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

第八章常微分方程的初值問題常微分方程（ODE）：只含有未知函數(shù)的導(dǎo)數(shù)的方程。ODE的階數(shù)：最高階導(dǎo)數(shù)的階數(shù)ODE問題初值問題：邊值問題：已知初始點(diǎn)處的條件已知初始點(diǎn)和末端點(diǎn)處的條件本章只討論初值問題一階ODE問題的形式1、如果f與y無關(guān)，則計(jì)算變?yōu)榈?章（數(shù)值積分）中討論的任一種直接積分方法應(yīng)用

初始條件2、如果f是y的函數(shù)，積分過程將不同于前者。若f是y的線性函數(shù)，如：f=ay+b其中a,b是常數(shù)或是t的函數(shù)，此時(shí)原方程稱為線性O(shè)DE若f不是線性函數(shù)，方程就稱為非線性O(shè)DE。一、求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à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例：求微分方程在初值條件下的特解，并畫出解函數(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)例：注4：解微分方程組時(shí)，如果所給的輸出個(gè)數(shù)與方程個(gè)數(shù)相同，則方程組的解按詞典順序輸出；如果只給一個(gè)輸出，則輸出的是一個(gè)包含解的結(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是一個(gè)結(jié)構(gòu)類型的數(shù)據(jù)r.x%查看解函數(shù)x(t)r.y%查看解函數(shù)y(t)dsolve的輸出個(gè)數(shù)只能為一個(gè)或與方程個(gè)數(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的一種簡(jiǎn)便方法。盡管精度不高，但由于簡(jiǎn)單，特別適用于快速編程求解。它有三種格式：二、用Euler方法求數(shù)值解（a）向前Euler法（b）改進(jìn)的Euler法（c）向后Euler法介紹這些方法是為了了解初值問題求解的基本思想。一階ODE對(duì)兩邊從x0到x積分得：(積分方程）1、向前Euler法推導(dǎo)1：設(shè)節(jié)點(diǎn)為向前Euler法：用向前差分公式代替導(dǎo)數(shù)：此處，y(xn)表示xn處的理論解，yn表示y(xn)的近似解一階ODE對(duì)兩邊從xn

到xn+1積分得：推導(dǎo)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法非常簡(jiǎn)單，但有兩種誤差效應(yīng)：第一種是截?cái)嗾`差，第二種誤差源于計(jì)算的不穩(wěn)定性。當(dāng)方程的時(shí)間系數(shù)為負(fù)值，h又不是很小時(shí)，第二種誤差就會(huì)增加。帶有負(fù)時(shí)間系數(shù)的常見方程為

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

到xn+1積分得：2、改進(jìn)的Euler法yn近似代替y(xn)用梯形代替右邊的積分梯形法從n=0開始計(jì)算，每步都要求解一個(gè)關(guān)于yn+1的方程（一般是一個(gè)非線性方程），可用如下的迭代法計(jì)算：利用此算法，可得：利用（為允許誤差）來控制是否繼續(xù)迭代迭代法太麻煩，實(shí)際上，當(dāng)h取得很小時(shí)，只讓上式中的第二式迭代一次就可以，即改進(jìn)的Euler法（也叫歐拉預(yù)估—校正法）預(yù)估算式校正算式改進(jìn)的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改進(jìn)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ù)，則與改進(jìn)的Euler算法一樣，在每一步中需要采用迭代法。該方法有兩個(gè)優(yōu)點(diǎn)：

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

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