MATLAB解微分方程

1、1MATLABMATLABODEODE初值問題的數(shù)值解初值問題的數(shù)值解PDEPDE問題的數(shù)值解問題的數(shù)值解2問題提出問題提出 倒葫蘆形狀容器壁上的刻度問題倒葫蘆形狀容器壁上的刻度問題. .對于如圖所示圓對于如圖所示圓柱形狀容器壁上的容積刻度柱形狀容器壁上的容積刻度, ,可以利用圓柱體體積公式可以利用圓柱體體積公式HDV22其中直徑其中直徑D D為常數(shù)為常數(shù). . 對于幾何形狀不是規(guī)則的容器對于幾何形狀不是規(guī)則的容器, ,比如倒葫蘆形狀容比如倒葫蘆形狀容器壁上如何標出刻度呢器壁上如何標出刻度呢? ?下表是經(jīng)過測量得到部分下表是經(jīng)過測量得到部分容器高度與直徑的關(guān)系容器高度與直徑的關(guān)系. .3x1o

2、根據(jù)上表的數(shù)據(jù)根據(jù)上表的數(shù)據(jù), ,可以擬合出倒可以擬合出倒葫蘆形狀容器的圖葫蘆形狀容器的圖, ,建立如圖所建立如圖所示的坐標軸示的坐標軸, ,問題為如何根據(jù)任問題為如何根據(jù)任意高度意高度x x標出容器體積標出容器體積V V的刻度的刻度, ,由微元思想分析可知由微元思想分析可知H 0 0.2 0.4 0.6 0.8 1.0D 0.04 0.11 0.26 0.56 1.04 1.17dxDdV241其中其中x x表示高度表示高度, ,直徑直徑D D是高度是高度x x的函數(shù)的函數(shù), ,記為記為D(xD(x).).x1o40)0()(412VxDdxdV只要求解上述方程只要求解上述方程, ,就可求出

3、體積就可求出體積V V與高度與高度x x之間的之間的函數(shù)關(guān)系函數(shù)關(guān)系, ,從而可標出容器壁上容積的刻度從而可標出容器壁上容積的刻度, ,但問題但問題是是函數(shù)函數(shù)D(xD(x) )無解析表達式無解析表達式, ,無法求出其解析解無法求出其解析解. . 因此因此, , 得到如下微分方程初值問題得到如下微分方程初值問題5包含自變量、未知函數(shù)及未知函數(shù)的導數(shù)或微分的包含自變量、未知函數(shù)及未知函數(shù)的導數(shù)或微分的方程稱為方程稱為微分方程微分方程。在微分方程中在微分方程中, , 自變量的個數(shù)只有一個自變量的個數(shù)只有一個, , 稱為稱為常微常微分方程。分方程。自變量的個數(shù)為兩個或兩個以上的微分方程叫自變量的個數(shù)

4、為兩個或兩個以上的微分方程叫偏微偏微分方程。分方程。微分方程中出現(xiàn)的未知函數(shù)最高階導數(shù)的階數(shù)稱為微分方程中出現(xiàn)的未知函數(shù)最高階導數(shù)的階數(shù)稱為微分方程的微分方程的階數(shù)階數(shù)。微分方程分類微分方程分類6常微分方程常微分方程：0)(),() 1 (yaybxayxfy )(,)(),()2(0ayyaybxayyxfy nybyyaybxayyxfy)(,)(),()3(02002212210012111)(),()(),()4(yxyyyxfyyxyyyxfy7但能求但能求解析解解析解的常微分方程是有限的，大多數(shù)的常的常微分方程是有限的，大多數(shù)的常微分方程是給不出解析解的微分方程是給不出解析解的.

5、. 22yxy 這個一階微分方程就不能用初等函數(shù)及其積這個一階微分方程就不能用初等函數(shù)及其積分來表達它的解。分來表達它的解。 例例 例例 1)0(yyy,xey 的解的解的值仍需插值方法來計算的值仍需插值方法來計算. .xe8222204,xxtdyxydxyee dt例 ：這是微積分的發(fā)明者之一Leibniz在1686年曾經(jīng)讓當時數(shù)學界人士求解的一階微分方程式，吸引了許多數(shù)學家的注意，大約經(jīng)過150年的探索，到1838年。劉維爾（Liouville)在理論上證明了這個微分方程不能用初等積分法求解，得借助于數(shù)值方法y =1-2xy例5：y(0)=0的解為為計算具體函數(shù)值，還需要數(shù)值積分方法。如

6、果需要許多點處的y值，則計算工作量較大?？梢?，用求解析的方法來找微分方程的數(shù)值解有時是不適宜的，必須研究微分方程的數(shù)值解法來求出其數(shù)值解9事實上，從事實上，從實際問題實際問題當中抽象出來的微分方程，通當中抽象出來的微分方程，通常主要依靠常主要依靠數(shù)值解法數(shù)值解法來解決。來解決。可以證明可以證明: : 如果函數(shù)在帶形區(qū)域如果函數(shù)在帶形區(qū)域 R=axb,R=axb,-y y內(nèi)連續(xù)，且關(guān)于內(nèi)連續(xù)，且關(guān)于y y滿足李普希茲滿足李普希茲(Lipschitz(Lipschitz) )條件，即存在常數(shù)條件，即存在常數(shù)L(L(它與它與x,yx,y無關(guān)無關(guān)) )使使 2121),(),(yyLyxfyxf對對R

7、 R內(nèi)任意兩個內(nèi)任意兩個 都成立都成立, ,則方程的解則方程的解 在在 a, ba, b 上存在且唯一。上存在且唯一。 21, yy)(xyy 00)(),(yxyyxfy在區(qū)間在區(qū)間a x ba x b上的數(shù)值解法上的數(shù)值解法。 主要討論一階常微分方程初值問題主要討論一階常微分方程初值問題 10微分方程數(shù)值方法的基本思想微分方程數(shù)值方法的基本思想 對常微分方程初值問題的數(shù)值解法，就是要算對常微分方程初值問題的數(shù)值解法，就是要算出精確解出精確解y(x)y(x)在區(qū)間在區(qū)間 a,ba,b 上的一系列離散節(jié)點上的一系列離散節(jié)點 處的函數(shù)值處的函數(shù)值bxxxxann110的近似 值)(,),(),(

8、10nxyxyxynyyy,10iixxh1niihxxi, 2 , 1 ,0相鄰兩個節(jié)點的間距相鄰兩個節(jié)點的間距 稱為步長，步稱為步長，步長可以相等，也可以不等。假定長可以相等，也可以不等。假定h h為定數(shù)，稱為定步長，為定數(shù)，稱為定步長，這時節(jié)點可表示為這時節(jié)點可表示為11數(shù)值解法需要把連續(xù)性的問題加以離散化，從數(shù)值解法需要把連續(xù)性的問題加以離散化，從而求出離散節(jié)點的數(shù)值解。而求出離散節(jié)點的數(shù)值解。離散化離散化 對常微分方程數(shù)值解法的基本出發(fā)點就是對常微分方程數(shù)值解法的基本出發(fā)點就是離散化離散化。其數(shù)值解法的基本特點：采用其數(shù)值解法的基本特點：采用“步進式步進式”，即求解過程順著節(jié)點排列的

9、次序一步一步地向前推進，即求解過程順著節(jié)點排列的次序一步一步地向前推進，12描述這類算法，要求給出用已知信息描述這類算法，要求給出用已知信息 計算計算 的遞推公式。的遞推公式。建立這類遞推公式的基本方法是在這些節(jié)點上用建立這類遞推公式的基本方法是在這些節(jié)點上用數(shù)數(shù)值積分、數(shù)值微分、泰勒展開值積分、數(shù)值微分、泰勒展開等離散化方法，對初等離散化方法，對初值問題值問題021,yyyyiii1iy00)(),(yxyyxfyy中的導數(shù)中的導數(shù) 進行不同的進行不同的離散化處理離散化處理。13數(shù)值解和精確解數(shù)值解和精確解用數(shù)值方法求解初值問題，不是求出它的解析解或其近似解析式，而是給出它的解在某些離散節(jié)點

10、上的近似值。用y(x)表示問題的準確解y(x0), y(x1),y(xN) 表示解y(x)在節(jié)點x0, x1, xN處的準確值y0,y1,y N表示數(shù)值解，即問題的解y(x) 在相應節(jié)點處的近似值。14單步法和多步法單步法和多步法單步法：在計算yi+1 時只利用y i多步法：在計算yi+1 時不僅利用y i , 還要利用 yi1, yi2,k步法：在計算yi+1 時要用到y(tǒng)i,yi1,yik+1顯式格式可寫成:yk+1=yk+hf(xk，yk;h）隱式格式:yk+1=yk+hf（xk，yk,yk+1;h）它每步求解yk+1需要解一個隱式方程。15歐拉（歐拉（EulerEuler）方法）方法在x

11、= x0 處，用差商代替導數(shù)：hxyxyxxxyxyxy)()()()()(0101010由由00000)(),()(yxyyxfxy10001),()(yyxhfyxy得得16同理，在同理，在x= xx= xn n 處，用差商代替導數(shù)：處，用差商代替導數(shù)：hxyxyxxxyxyxynnnnnnn)()()()()(111),()()(1nnnnyxhfxyxy),()(nnnyxfxy由由得得若記若記11)(,)(nnnnyxyyxy則上式可記為則上式可記為),(1nnnnyxhfyy17例例: : 用用EulerEuler方法求解常微分方程初值問題方法求解常微分方程初值問題 yyxyxy2

12、03002()( ).并將數(shù)值解和該問題的解析解比較。并將數(shù)值解和該問題的解析解比較。21)(xxxy解析解：解：解：EulerEuler方法的具體格式：方法的具體格式：)2(21nnnnnyxyhyy18h=0.2;y(1)=0.2;x=0.2:h:3;h=0.2;y(1)=0.2;x=0.2:h:3;for n=1:14for n=1:14 xn=x(n);yn xn=x(n);yn=y(n);=y(n); y(n+1)=yn+h y(n+1)=yn+h* *(yn/xn-2(yn/xn-2* *ynyn* *ynyn););endendx0=0.2:h:3;y0=x0./(1+x0.2)

13、;x0=0.2:h:3;y0=x0./(1+x0.2);plot(x0,y0,x,y,x,y,o)plot(x0,y0,x,y,x,y,o)程序?qū)崿F(xiàn)程序?qū)崿F(xiàn)19xn y(xn) yn yn-y(xn)0.00000.20.19230.20000.00770.40.34480.38400.03920.60.44120.51700.07580.80.48780.58240.09461.00.50000.59240.09241.20.49180.57050.07871.40.47300.53540.0624h h=0.2, =0.2, x xn n= =nhnh,(,(n n=0,1,2=0,1,2

14、,15), ,15), f f( (x x, ,y y)=)=y y/ /x x 2 2y y2 2 計算中取計算中取f f(0,0)=1. (0,0)=1. 計算結(jié)果如下：計算結(jié)果如下：20 xn y(xn) yn yn-y(xn)1.60.44940.49720.04781.80.42450.46050.03592.00.40000.42680.02682.20.37670.39660.01992.40.35500.36980.01472.60.33510.34590.01082.80.31670.32460.00793.00.30000.30570.0057由表中數(shù)據(jù)可以看到，微分方程初

15、值問題的數(shù)值解和解由表中數(shù)據(jù)可以看到，微分方程初值問題的數(shù)值解和解析解的誤差一般在小數(shù)點后析解的誤差一般在小數(shù)點后第二位或第三位小數(shù)第二位或第三位小數(shù)上，說上，說明明EulerEuler方法的精度是比較差的。方法的精度是比較差的。21二階二階Runge-KuttaRunge-Kutta方法方法),(),()(12122122111hKyhxfKyxfKKcKchyynnnnnn),(),()(2122111nnnnnnnnhfyhxfcyxfchyxyT22即即021021012122221cccc212112122221cccc23例例 蛇形曲線蛇形曲線的初值問題的初值問題令令f(x,y)=

16、y/x 2yf(x,y)=y/x 2y2 2, , 取取 f(0,0)=1, f(0,0)=1, h=0.2,xh=0.2,xn n=hn=hn , ( n = 1,2, , ( n = 1,2,15),15)2 2階龍格階龍格- -庫塔公式計算格式：庫塔公式計算格式： k k1 1=y=yn n/x/xn n 2y 2yn n 2 2, , k k2 2 = (y = (yn n+hk+hk1 1)/(x)/(xn n+h) +h) 2(y2(yn n+hk+hk1 1) )2 2 y yn+1n+1=y=yn n + 0.5h k + 0.5h k1 1 + k + k2 2 22( 03

17、 )(0)0.yyyxxy 2400.511.522.5x0=0;y0=0;h=.2;x=.2:h:3;k1=1;k2=(y0+h*k1)/x(1)-2*(y0+h*k1)2;y(1)=y0+.5*h*(k1+k2);for n=1:14 k1=y(n)/x(n)-2*y(n)2; k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)2; y(n+1)=y(n)+0.5*h*(k1+k2);endy1=x./(1+x.2);plot(x,y,o,x,y1)25),()2,2()2,2(),()22(6342312143211hKyhxfK

18、KhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn四階四階Runge-KuttaRunge-Kutta方法方法26)(txx )()()(txtxty)()()()()(12tytytxtxty function ydot = harmonic(t,y) ydot=y(2); -y(1) y=inline(0 1;-1 0*y,t,y);System of Equations2733)(/ )()()(/ )()(trtvtvtrtutu 22)()()(tvtutr)()()()()(tvtutvtuty33)(/ )()(/ )()()()(trtvtrtutvtuty

19、 function ydot =twobody(t,y) r=sqrt(y(1)2+y(2)2); ydot=y(3);y(4);-y(1)/r3;-y(2)/r3; Two Body Problem28)()(),(),(ccccyyJttytfytf),(),(ccccytyfJyttf),.,(),.,(),.,()()()(1121121nnnnnyytfyytfyytftytytydtdnnnnnnyfyfyfyfyfyfyfyfyfJ212221212111Jyy Linearized Differential Equations29J的特征值是kkki)(diagk1VVJyVx

20、 kkkxx)()()(cttktxetxck0k0k0k解增長解衰減解振蕩30 基于龍格庫塔法，基于龍格庫塔法， MATLABMATLAB求常微分方程數(shù)值解求常微分方程數(shù)值解的函數(shù)的函數(shù)，一般調(diào)用格式為：一般調(diào)用格式為： t,yt,y=ode23(fname,tspan,y0)=ode23(fname,tspan,y0) t,y t,y=ode45(fname,tspan,y0)=ode45(fname,tspan,y0)其中其中fnamefname是定義是定義f(t,yf(t,y) )的函數(shù)文件名的函數(shù)文件名，該函數(shù)文，該函數(shù)文件必須返回一個列向量。件必須返回一個列向量。tspantspa

21、n形式為形式為t0,tf,t0,tf,表表示求解區(qū)間。示求解區(qū)間。y0y0是初始狀態(tài)列向量。是初始狀態(tài)列向量。t t和和y y分別給分別給出時間向量和相應的狀態(tài)向量。出時間向量和相應的狀態(tài)向量。MATLABMATLAB求常微分方程數(shù)值解的函數(shù)求常微分方程數(shù)值解的函數(shù)31 ode23: Bogacki,Shampine(1989)和Shampine(1994), ”23”表示用兩算法：一個2階，一個3階)9865(72),()432(9)43,43()2,2(),(432111143211123121ssssheytfsssshyyhtthsyhtfsshyhtfsytfsnnnnnnnnnnn

22、nnBogacki, P. and LF Shampine, A 3(2) pair of Runge-Kutta formulas, Appl.Math. Letters, Vol. 2, 1989, pp 1-9. BS23 algorithm32 F=inline(y(2);-y(1),t,y) ode23(F,0 2*pi,1;0) opts=odeset(reltol,1.e-4,abstol,1.e-6,outputfcn)Examples33 ode23(twobody,0 2*pi,1;0;0;1);01234567-1.5-1-0.500.511.5Examples34 y0

23、=1;0;0;3; ode23(twobody,0 2*pi,y0); 01234567-202468101214161835 y0=1;0;0;3; t,y=ode23(twobody,0 2*pi,y0); plot(y(:,1),y(:,2); axis equal-1.5-1-0.500.5102468101214161836 y0=1;0;0;3; t,y=ode23(twobody,0 2*pi,y0); plot(y(:,1) plot(y(:,2)37 A problem is stiff if the solution being sought is varying slow

24、ly, but there are nearby solutions that very rapidly, so the numerical method must take small steps to obtain satisfactory results.A model of flame propagation（火焰燃燒）:/20)0(32tyyyy y是球的半徑，y2和y3與球的表面積和體積有關(guān)想一下，點燃一根火柴，光球迅速增長，到達一個關(guān)鍵的大小，然后維持它的大?。ㄓ捎谶M入球內(nèi)氧氣和消耗的氧氣平衡）Stiff Problem(剛性問題)380102030405060708090100

25、00.811.21.4eta=0.02;sym y; F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);39eta=0.00002;sym y; F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);40012345678910 x 10400.811.21.4ode23seta=0.00002; ode23s(inline(y2-y3,t,y),0 2/eta,eta);41例例 蛇形曲線的常微分方程初值問題蛇形曲線的常微分方程初值問題 22211yxy 0)0( yMATLABMATLAB數(shù)值求

26、解命令數(shù)值求解命令F=inline(1./(1+x.2)-2F=inline(1./(1+x.2)-2* *y.2)y.2)；ode23(F,0,6,0)ode23(F,0,6,0)輸出結(jié)果為圖形輸出結(jié)果為圖形 42T,y=ode23(f,0,6,0)T,y=ode23(f,0,6,0)將將得到自變量和函數(shù)的離散數(shù)得到自變量和函數(shù)的離散數(shù)據(jù)據(jù) T =0 0.0001 0.0005 0.0025 0.0125 0.0496 0.1085 0.1863 0.2837 0.4091 0.5991 0.8513 1.0567 1.2680 1.5110 1.8050 2.1788 2.6842 3.2

27、842 3.8842 4.4842 5.0842 5.6842 6.0000y =0 0.0001 0.0005 0.0025 0.0125 0.0495 0.1073 0.1800 0.2626 0.3505 0.4411 0.4944 0.5000 0.4868 0.4607 0.4242 0.3793 0.3270 0.2783 0.2411 0.2122 0.1891 0.1705 0.162043例例 洛倫茲模型洛倫茲模型由如下常微分方程組描述由如下常微分方程組描述 zyxydtdzzydtdyyzxdtdx )(取取 = =8/38/3， = =1010， = =2828。初值初值

28、：x x(0)=0(0)=0，y y(0)=0(0)=0，z z(0)=0.01(0)=0.01。利用利用MATLABMATLAB求解常微分方程數(shù)值解命令計算出求解常微分方程數(shù)值解命令計算出t t00，8080內(nèi)內(nèi)，三個未知函數(shù)的數(shù)據(jù)值，并繪出相空間在三個未知函數(shù)的數(shù)據(jù)值，并繪出相空間在Y-XY-X平面的投影曲線平面的投影曲線 44氣象學家Lorenz提出一篇論文，名叫一只蝴蝶拍一下翅膀會不會在Taxas州引起龍卷風？論述某系統(tǒng)如果初期條件差一點點，結(jié)果會很不穩(wěn)定，他把這種現(xiàn)象戲稱做蝴蝶效應蝴蝶效應。Lorenz為何要寫這篇論文呢？這故事發(fā)生在1961年的某個冬天，他如往常一般在辦公室操作氣象

29、電腦。平時，他只需要將溫度、濕度、壓力等氣象數(shù)據(jù)輸入，電腦就會依據(jù)三個內(nèi)建的微分方程式，計算出下一刻可能的氣象數(shù)據(jù)，因此模擬出氣象變化圖。這一天，Lorenz想更進一步了解某段紀錄的後續(xù)變化，他把某時刻的氣象數(shù)據(jù)重新輸入電腦，讓電腦計算出更多的後續(xù)結(jié)果。當時，電腦處理數(shù)據(jù)資料的數(shù)度不快，在結(jié)果出來之前，足夠他喝杯咖啡并和友人閑聊一陣。在一小時後，結(jié)果出來了，不過令他目瞪口呆。結(jié)果和原資訊兩相比較，初期數(shù)據(jù)還差不多，越到後期，數(shù)據(jù)差異就越大了，就像是不同的兩筆資訊。而問題并不出在電腦，問題是他輸入的數(shù)據(jù)差了0.000127，而這些微的差異卻造成天壤之別。所以長期的準確預測天氣是不可能的。45天氣

30、預報的準確性：http:/.tw/senior/seniorteach/stfastnews/stnfastnews/stnfastnews4/stnfastnews4_2.htmLorenz現(xiàn)象的數(shù)學： /news_detail.php?news_id=225分形藝術(shù)電子版：http:/ y y1 1，y y2 2，y y3 3 = = x x，y y，z z ，創(chuàng)建創(chuàng)建MATLABMATLAB函函數(shù)文件如下數(shù)文件如下 function z=flo(t,yfunction z=flo(t,y) )z(1,:)=-8z(1,:)=-8* *y(1)/

31、3+y(2).y(1)/3+y(2).* *y(3);y(3);z(2,:)=-10z(2,:)=-10* *(y(2)-y(3);(y(2)-y(3);z(3,:)=-y(1).z(3,:)=-y(1).* *y(2)+28y(2)+28* *y(2)-y(3);y(2)-y(3);用用MATLABMATLAB命令求解并繪出命令求解并繪出Y-XY-X平面的投影圖平面的投影圖 y0=0;0;0.01;y0=0;0;0.01;x,y=ode23(flo,0, 80,y0)x,y=ode23(flo,0, 80,y0)plot(y(:,2),y(:,1) plot(y(:,2),y(:,1) 47

32、-20-15-10-50510152001020304050plot(y(:,3),y(:,1) plot(y(:,3),y(:,1) 48plot3(y(:,1),y(:,2),y(:,3)49 y0=30;0;-40; plot(y(:,iy0=30;0;-40; plot(y(:,i)505152非剛性系統(tǒng)： ode45(Runge-Kutta45) ode23(Runge-Kuatta23) ode113(Adams-Bashforth-Moulton PECE)多步方法剛性系統(tǒng)： ode15s(Gear方法) ,多步方法 ode23s(二階modified Rosenbroack f

33、ormula)，單步 ode23t(trapezoidal rule),solve DAEs ode23tb(TR-BDF2) low order methodMatlabs ODE Solvers532222yxLaplacian 算子： Poisson方程（elliptic）：fu Laplacian 算子的特征值問題:fuuHeat equation(parabolic):utuWave equation(hyperbolic):utu22PDE Model54五點離散22),(),(2),( ),(),(2),(),(hhyxuyxuhyxuhyhxuyxuyhxuyxuh2)(4)(

34、)()()()(hPuSuEuWuNuPuh0)(Puh)()(PfPuh)()(PuPukkh Poisson方程離散： 特征值問題:Finete Difference Methods55),(),(),(yxutyxutyxuh),(),(),(yxutyxutyxuh),(),(),(2),(2yxutyxutyxutyxuh),(),(),(2),(2yxutyxutyxutyxuh熱方程：波動方程：56-4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 -4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 -4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -4