Matlab解微分方程(ODE+PDE)

1、常微分方程:1 ODE解算器簡介(ode*)2微分方程轉(zhuǎn)換3剛性/非剛性問題(Stiff/Nonstiff)4隱式微分方程(IDE)5微分代數(shù)方程(DAE)6延遲微分方程(DDE)7邊值問題(BVP)偏微分方程(PDEs)Matlab解法偏微分方程：1 一般偏微分方程組(PDEs)的命令行求解2特殊偏微分方程(PDEs)的PDEtool求解3陸君安偏微分方程的 MATLAB解法先來認識下常微分方程(ODE)初值問題解算器(solver)T,Y,TE,YE,IE = odesolver(odefun,tspan,yO,optionS sxint = deval(sol,xint)Matlab中提

2、供了以下解算器:解曽樂udEoivtr |趕44?* 5非m hrranztiffln：T- 4hi-憐式” Afta *ihiirif dr 2 NJ L. .tjLT K:=_耳.耳時T.豈x比ihiPtns.+ *.串K電十鬥亍卄韭用煒更曲拼狷拼 心iSJTAuu重"麒lUJT展悝可"RJT廿匕上”慎現(xiàn) r+hrfii王古呼護離.占就茯用叫権常f. r九7-:*燉匕.確半囲寸馬馬H苕埼再+ 計/叭ME In 'i： ± 7-障佇ft ft+±itfS住訊翳的聚鴉翠肝僉式，兇星沖田1*好什策 tW+ RClttHITmCft歯千言彈甘篇+ d

3、tWltJRtttS悔些.r|t. <r.v£如14輸入?yún)?shù):Matlabodefun：微分方程的Matlab語言描述函數(shù)，必須是函數(shù)句柄或者字符串，必須寫成規(guī)范格式 (也就是一階顯示微分方程組)，這個具體在后面講解tspan=tO tf或者tO,t1,tf微分變量的范圍，兩者都是根據(jù)tO和tf的值自動選擇步長，只是前者返回所有計算點的微分值，而后者只返回指定的點的微分值，一定要注意對于后者tspan 必須嚴格單調(diào) ，還有就是兩者數(shù)據(jù)存儲時使用的內(nèi)存不同(明顯前者多 )，其它沒有任何本質(zhì)的區(qū)別y0=y(0),y(0),y '：微分方程初值，依次輸入所有狀態(tài)變量的初值，什

4、么是狀態(tài)變量在后面有介紹opti ons：微分優(yōu)化參數(shù)，是一個結(jié)構(gòu)體，使用odeset可以設(shè)置具體參數(shù)，詳細內(nèi)容查看幫助輸出參數(shù)：T：時間列向量，也就是 ode*計算微分方程的值的點Y :二維數(shù)組，第i列表示第i個狀態(tài)變量的值，行數(shù)與T 一致在求解ODE時，我們還會用到deval()函數(shù)，deval的作用就是通過結(jié)構(gòu)體solution計算t對應(yīng) x 值，和 polyval 之類的很相似！參數(shù)格式如下：sol:就是上次調(diào)用 ode*函數(shù)得道的結(jié)構(gòu)體解xint:需要計算的點，可以是標(biāo)量或者向量，但是必須在tspan范圍內(nèi)該函數(shù)的好處就是如果我想知道 t=tO 時的 y 值，不需要重新使用 ode

5、計算，而直接使用上 次計算的得道 solution 就可以 教程 微分方程轉(zhuǎn)換為一階顯示微分方程組方法好，上面我們把 Matlab中的常微分方程(ODE)的解算器講解的差不多了，下面我們就具體 開始介紹如何使用上面的知識吧！現(xiàn)實總是殘酷的，要得到就必須先付出，不可能所有的 ODE-拿來就可以直接使用，因此, 在使用ODEW算器之前，我們需要做的第一步，也是最重要的一步， 借助狀態(tài)變量將微分方程組化成Matlab可接受的標(biāo)準形式（一階顯示常微分方程）!如果ODEs由一個或多個高階微分方程給出，則我們應(yīng)先將其變換成一階顯式常微分方程組！ 下面我們以兩個高階微分方程構(gòu)成的ODEs為例介紹如何將之變換

6、成一個一階顯式常微分方程組。stepl.將微分方程的最高階變量移到等式的左邊，其他移到右邊，并按階次從低到高排列，假如說兩個高階微分方程最后能夠顯式的表達成如下所示：/ ,滬7 J J；$7）我們說過現(xiàn)實總是殘酷的，有時方程偏偏是隱式的，沒法寫成上面的樣子， 不用擔(dān)心Matlab早就為我們想到了，這個留在后面的隱式微分方程數(shù)值求解中再詳細講解！step2.為每一階微分式選擇 狀態(tài)變量，最高階除外xm+i -y* xn+ y >> 忌從上面的變換，我們注意到，ODEs中所有因變量的最高階次之和就是需要的狀態(tài)變量的個數(shù)，最高階的微分式（比如上面的x（m）和y（n）不需要給它一個狀態(tài)變量

7、step3.根據(jù)上面選用的狀態(tài)變量，寫出所有狀態(tài)變量的一階微分的表達式X 賊+JS注意到，最高階次的微分式的表達式直接就是stepl中的微分方程好，到此為止，一階顯式常微分方程組，變化順利結(jié)束，接下來就是Matlab編程了。請記住上面的變化很重要，Matla中所有微分方程的求解都需要上面的變換。F面通過一個實例演示 ODEs的轉(zhuǎn)換和求解己卸Apoll衛(wèi)星的運動軌跡（兀刃滿足下面的槪分方程組= 2y'+ x-y " = -2x'-i-y譏并+旳乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高階次的微分式的表達式直接就是stepl中的微分方程乳4 =

8、 = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高階次的微分式的表達式直接就是stepl中的微分方程其中，q = Jd + nF+yF = y/(x-y +/,出二打 82 45x(0) = 1.2,x'(0) = 0JXO) = 0 j"(0) = -l 04935751【解】 真是萬幸，該 ODEs已經(jīng)幫我們完成了 stepl，我們只需要完成 step2和step3了（1）選擇一組狀態(tài)變量乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高階次的微分式的表達式直接就是stepl中的微分方程乳4 = = 2兀2 + 兀3 口 JTj / q /

9、 乜注意到，最高階次的微分式的表達式直接就是stepl中的微分方程（2）寫出所有狀態(tài)變量的一階微分表達式乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高階次的微分式的表達式直接就是stepl中的微分方程乳4 = = 2兀2 + 兀3 口 JTj / q / 乜注意到，最高階次的微分式的表達式直接就是stepl中的微分方程2 = 2百+龍】一 “（可+“）r三”（兀一“）/勺乳4 = = 2兀2 + 兀3 口 JTj / q / 乜有了數(shù)學(xué)模型描述，則可以立即寫出相應(yīng)的Matlab代碼了（這里我們優(yōu)先選擇 ode45）1. x0=1.2;0;0;-1.04935751;%x

10、0(i)對應(yīng)與 xi的初值2. options=odeset('reltol',1e-8);% 該命令的另一種寫法是options=odeset;options.reltol=1e-8;3. tic4. t,y=ode45(appollo,0,20,x0,options);%t 是時間點，y 的第 i 列對應(yīng) xi 的值，t 和 y 的行數(shù)相同5. toc6. plot(y(:,1),y(:,3)%繪制x1和x3,也就是x和y的圖形7. title('Appollo衛(wèi)星運動軌跡')8. xlabel( X)9. ylabel('Y')10.10.

11、 function dx=appollo(t,x)11. mu=1/82.45;12. mustar=1-mu;13. r1= sqrt(x(1)+mu)A2+x(3)A2);14. r2=sqrt(x(1)-mustar)A2+x(3)A2);15. dx=x(2)16. 2*x(4)+x(1)-mustar*(x(1)+mu)/r1A3-mu*(x(1)-mustar)/r2A317. x(4)18. -2*x(2)+x(3)-mustar*x(3)/r1A3-mu*x(3)/r2A3;教程剛性/非剛性問題（Stiff/Nonstiff）的Matlab解法在工程實踐中，我們經(jīng)常遇到一些OD

12、Es其中某些解變換緩慢，另一些變化很快，且相差懸殊的微分方程，這就是所謂的剛性問題（Stiff），對于所有解的變化相當(dāng)我們則稱為非剛性問題（Nonstiff）。對于剛性問題一般不適合使用ode45這類函數(shù)求解。由于非剛性問題我們使用的多比較多，我們就不多說，下面主要講解下Stiff看一個例子，考慮下面的微分方程” =0 04（1-孫）-（1-必）” + 0.0001（1-）2=-10000 + 30000-其中乃（0）二片3）二1 OJOO【解】 題目中已經(jīng)幫我們完成了方程組轉(zhuǎn)換，下面我們就直接編程了（1）編寫微分方程函數(shù)1. odefun=(t,x)0.04*(1-x(1)-(1-x(2)*

13、x(1)+0.0001*(1-x(2)F22. -1e4*x(1)+3000*(1-x(2)F2;3. x0=0 1;4. tspan=0 100;太慢了)5. options=odeset('reltol',1e-6,'abstol',1e-8);(2) 對于這個剛性問題，我們先使用ode45函數(shù)試試(反正我們是沒有信心等待它,1. tic;t,y=ode45(odefun,tspan,x0,options);toc2. disp('ode45 計算的點數(shù)'num2str(length(t) 換用ode15s函數(shù)試試看看1. tic;t2,y2

14、=ode15s(odefun,tspan,x0,options);toc2. disp('ode15s 計算的點數(shù)'num2str(length(t2)(4)繪圖看看結(jié)果1. figure('numbertitle','off','name','Stiff ODEs Demo by Matlabsky')2. subplot(121)3. title('Using ode45')4. plot(t,y)5. subplot(122)6. plot(t2,y2)7. title('Using

15、ode15s')運行的結(jié)果如下，我們可以比較下ode45和ode15s的差距1. Elapsed time is 171.005688 seconds.2. ode45計算的點數(shù)3569813. Elapsed time is 0.397287 seconds.4. ode15s計算的點數(shù)1881L,!Srng odsX&1U&mg心他甜-0?0B-C.B0.7-0.?2OS-0.50.S卩t-C403-03*0?-0.2010.1J-r.D°(U1G08教程隱式微分方程(IDE)的Matlab解法上帝不會總是那么仁慈的，不是所有的ODEs都是可以直接顯式的表

16、達成下面的樣子比如，下面的隱式微分方程組才+沙6"l>l+3x>"+zyl-z = 5那該如何辦呢?在這里我們介紹三種解決方法，但是前兩者是技巧方法，沒有使用Mathworks專為IDE開發(fā)的ode15i函數(shù)：使用符號計算函數(shù) solve()求解出x''和y''的表達式，不就可以了嗎？那好，我們下面試試看看解(1)求解表達式1.2. >> d2x,d2y=solve('d2x+2*dy*x-2*d2y','d2x*dy+3*dx*d2y+x*dy-y-5','d2x',&

17、#39;d2y')%d2x 表示 x 的二階導(dǎo)數(shù)，其它同理3.3. d2x =5.4. -2*(3*dy*x*dx-5+dy*x-y)/(3*dx+2*dy)7.5. d2y =9.6. (2*dyA2*x+5-dy*x+y)/(3*dx+2*dy)也就是說原來的ODEs可以重新寫成亍=_2(3''-5+7-7)3忑'+2”y =也就是說stepl完成了，接下來的就比較簡單了(2).選取狀態(tài)變量可二兀抵二X：礙二珂二卩(3) .寫出狀態(tài)變量的一階微分表達式xj =嚴=_2張存1眄J+旺-可)(4) .寫Matlab程序進行數(shù)值求解1. myode=(t,x)x(

18、2)2. -2*(3*x(1)*x(2)*x(4)-5+x(4)-x(3)/(3*x(2)+2*x(4)3. x(4)4. (2*x(4)A2*x(1)+5-x(4)*x(1)+x(3)/(3*x(2)+2*x(4)5. t,y=ode45(myode,0,20,x0);但是solve()函數(shù)也不是萬能的，對有些完全隱式，solve()可能壓根沒法求解出它們的具 體表達式，那此時該怎么辦呢？天無絕人之路，車到山前必有路，下面我們分析下數(shù)值解 法的本質(zhì)。其實Matlab要求我們先將ODEs轉(zhuǎn)化為顯式的一階微分方程組，就是為了便于計算狀態(tài)變量一階微分函數(shù)值(注意只需要它的值)！雖然現(xiàn)在我們沒有解析

19、的得到各個狀態(tài)變量一階 微分的表達式，但是我們只要求出每個點的對應(yīng)的所有狀態(tài)變量一階微分值即可。這個思想就是在 Matlab沒有提供ode15i()函數(shù)之前，大家唯一能使用的方法，下面我們以個例子說明下fsiii。)十 =-2 砂Ixx'*y % cosfj;") = 3yx'其中xO=XQO;l對于這樣的ODEs難道你還想將它化為那個一階顯式微分方程組嗎，別費勁了 ！【解】(1) 好，我們同樣的需要一組狀態(tài)變量西二兀 可二瑞 兀二兒兀=yl(2) .寫出狀態(tài)變量的一階微分。 只不過此時的x''和y''不再是顯式的，我們使用fsolve

20、進行數(shù)值求解花'sin() + (x4 )3 = 一細碼 + 碼花 H珂；vj 筍 + cos(x4 *) =編寫Matlab代碼1. t,y=ode45(odefun,03,1 0 0 1');2. plot(t,y)3. Iegend('x','x'",'y',y”)4.4. function dy=odefun(t,x)5. fun=(dxy)dxy(1)*sin(x(4)+dxy(2F2+2*x(1)*x(3)-x(1)*dxy(1)*x(4)6. x(1)*dxy(1)*dxy(2)+cos(dxy(2)-3*

21、x(3)*x(2);7. options=optimset('display','off');y=fsolve(fun,x(1,3),options);%使用fsolve求解出x"和 y"10.dy=x(2);y(1);x(4);y(2);% 狀態(tài)變量一階微分值復(fù)制代碼方法二也有它的缺點，每一個點的x“和y''都需要獨立計算， 假如說數(shù)據(jù)量很大的話，那個叫難熬呀，并且有時可能是由于方程或者參數(shù)配置問題fsolve還有會出現(xiàn)罷工現(xiàn)象！于是MathWorks為我們準備了 ode15i函數(shù)，不過這個 ode15i有點不好用，沒有 o

22、de45等 那么直接。ode15i()調(diào)用格式如下yOmod, yp0mod,resnrm = decic(odefun,tO,yO,fixed_yO,ypO,fixed_ypO,options.)T,Y,TE,YE,IE = ode15i(odefu n,tspa n, y0mod,yp0mod,optio ns.)隱式ODEs不同于一般的顯式 ODEs ode15i需要同時給出狀態(tài)變量 及其一階導(dǎo)數(shù)的初值（xO,dxO），它們不能任意賦值，只能有n個是獨立的，其余的需要隱式方程求解（使用decic 函數(shù)），否則出現(xiàn)矛盾的初值條件。輸入?yún)?shù)：odefun :微分方程，注意此時的函數(shù)變量有t（

23、自變量），x（狀態(tài)變量），dx（狀態(tài)變量一階導(dǎo)數(shù)）等三個tO :自變量初值，必須與 tspan 相等y0 :狀態(tài)變量初值，題目給出幾個就寫幾個，沒給的自己猜測一個fixed_yO :與yO的尺寸一樣，表示對應(yīng)位置的初值是給定的還是猜測的，如果給定則1，否則0ypO:狀態(tài)變量一階導(dǎo)數(shù)初值，其它同yOfixed_ypO :同理fixed_yO，不過此時對應(yīng)的是 ypO 了options :其它選項，具體查看幫助tspan :自變量的范圍，其中 tspan（1）必須與tO 一致輸出參數(shù)：大部分同理ode45，這里就不具體說明了F面我們以實例說明，對于下面的IDE，求解x"£1口0

24、) + (”呎=-2xy-xxny' "j'+cosO-) = 3'其中zO = l;ftO;l解（1）選取狀態(tài)變量西二兀 叼二瑞 為二兒旺=y(2)寫出狀態(tài)變量的一階微分。只不過此時的x''和y''不再是顯式的Xj '5111(X4)4-尸=-2磯可十兀/J 為珂花久*)=鄉(xiāng)也也書寫odefun，注意此時多了一個 dx變量，就是x的一階導(dǎo)數(shù)1. odefun=(t,x,dx)dx(1)-x(2)2. dx(2)*si n(x(4)+dx(4F2+2*x(1)*x(3)-x(1)*dx(2)*x(4)3. dx(3)-x

25、(4)4. x(1)*dx(2)*dx(4)+cos(dx(4)-3*x(3)*x(2);(4)求解dx01. t0=0;%自變量初值2. %對于x0和dx0中的題目給出的初值，如實寫，沒有給出的任意寫3. x0=1 0 0 1'%本題初值x0的都給出了，所以必須是這個4. %初值x0中題目給出的，對應(yīng)位置填1,否則為05. fix_x0=ones(4,1);%本題中x0都給出了，故全為 16. dx0=0 0 1 1'%本題中初值dx0 一個都沒有給出，那么全部任意寫7. %初值dx0中題目給出的，對應(yīng)位置填 1，否則為08. fix_dx0=zeros(4,1);%本題中d

26、x0 一個沒有給出，故全部為 09. x02,dx02=decic(odefun,t0,x0,fix_x0,dx0,fix_dx0);(5) 求解微分方程1. solution=ode15i(odefun,0 2,x02,dx02)%x02 和 dx02 是由 decic 輸出參數(shù)雖然我們上面的分析是完全正確，但是好像Matlab就是罷工，郁悶死了其實對一些簡單的 ODEs壓根沒有必要通過 decic函數(shù)來求解未知的初值，我們可以直接人工算出x02和dx02，根據(jù)下面的式子花 011（駐）+ E 忙=一2誥內(nèi) + X2 *X4我們演示下，根據(jù)題目x0=x1,x2,x3,x4=1 0 0 1，由

27、第一個式子可以得出x1' =0,第三個式子有x3' =1，再聯(lián)立第二四兩個式子后有x2' =0，x4' =1，故dx0=x1 ' ,x2 ' ,x3 ' ,x4 ' =0 0 1 1看了下，隱式 ODEs數(shù)值求解的確麻煩的些，但是和前面的ode45沒有本質(zhì)的區(qū)別，decic函數(shù)就是為了求解那些隱式中的狀態(tài)變量的一階導(dǎo)數(shù)，與我們方法二中的很相似，只不過 這里直接調(diào)用罷了。其實不管隱式還是顯式 ODEs甚至DAE都是可以使用 ode15i求解，只不過此時將顯式當(dāng) 隱式處理，或者說將代數(shù)約束看成一個隱式罷了！這一條等到我們在介紹 DAE

28、的時候再驗證！當(dāng)然大部分人都不會傻到這么做，有現(xiàn)成的函數(shù)不用，何必自找麻煩呢！但是對于了解和 學(xué)習(xí)，我們可以試試。教程微分代數(shù)方程(DAE)的Matlab解法所謂微分代數(shù)方程，是指在微分方程中，某些變量滿足某些代數(shù)方程的約束。假設(shè)微分方 程的更一般形式可以寫成''''''前面所介紹的ODEs數(shù)值解法主要針對能夠轉(zhuǎn)換為一階常微分方程組的類型，故DAE就無法使用前面介紹的常微分方程解法直接求解，必須借助DAE的特殊解法。其實對于我們使用 Matlab求解DAE時，卻沒有太大的改變 只需增加一個 Mass參數(shù)即可。 描述f(t,x)的方法和普通微分方程

29、完全一致。、，、，、廠注意： ode15i 沒法設(shè)置 Mass 屬性，換句話說除了 ode15i 外其他 ode 計算器都可以求解DAEs問題1. 當(dāng) M(t,y) 非奇異的時候， 我們可以將微分方程等效的轉(zhuǎn)換為 y'=inv(M)*f(t,y) ，此時就是一個普通的ODE當(dāng)然我們可以將它當(dāng)成 DAEs處理)，對任意一個給定的初值條件都有唯一的解2. 當(dāng)m(t,y)奇異時，我們叫它為 DAEs(微分代數(shù)方程)，DAEs問題只有在同時提供狀態(tài)變量初值y0和狀態(tài)變量一階導(dǎo)數(shù)初值pyO，且滿足M(tO,yO)*ypO=f(tO,yO)時才有唯一解，假如不滿足上面的方程，DAEs解算器會將提供

30、的yO和pyO作為猜測初始值，并重新計算與提供初值最近的封閉初值3. 質(zhì)量矩陣可是一個常數(shù)矩陣 (稀疏矩陣 )，也可以是一個自定義函數(shù)的輸出。但是 ode23s只能求解Mass是常數(shù)的DAEs4. 對于Mass奇異的DAEs問題，特別是設(shè)置 MassSingular為yes時，只能ode15s和ode23t 解算器5. 對于DAE我們還有幾個參數(shù)需要介紹a. Mass :質(zhì)量矩陣，不說了，這個是DAE的關(guān)鍵，后面看例子就明白了b. MStateDependence ：質(zhì)量矩陣 M(t,y) 是否是 y 的函數(shù)，可以選擇 none|weak|strong ， none 表示 M與 y 無關(guān)，we

31、ak禾口 strong 者E表示與 y 相關(guān)c. MvPattern :注意這個必須是稀疏矩陣， S(i,j)=1 表示 M(t,y) 的第 i 行中任意元素都與 第 j 個狀態(tài)變量 yi 有關(guān)，否則為 Od. MassSingular :設(shè)置 Mass矩陣是否奇異，當(dāng)設(shè)置為yes時，只能使用 ode15s和ode23te. InitialSlope:狀態(tài)變量的一階導(dǎo)數(shù)初值yp0，和y0具有相同的size，當(dāng)使用ode15s和 ode23t 時，該屬性默認為 O面我們以實例說明，看下面的例子，求解該方程的數(shù)值解珂 =0-2s)十 叼珂 十0 2珂屯 ' = 2屯;t - 5x32r?

32、- 2才再+吃+旳=1其中匝(°) 0 & 旳(°)=忑(0) = 0.1解真是萬幸，選取狀態(tài)變量和求狀態(tài)變量的一階導(dǎo)數(shù)等，微分方程轉(zhuǎn)換工作，題目已經(jīng)幫我 們完成。DAE可是細心的網(wǎng)友會發(fā)現(xiàn)，最后一個方程不是微分方程而是一個代數(shù)方程(這就是為什叫的原因)，其實我們可以將它視為對三個狀態(tài)變量的約束。(1) 用矩陣形式表示出該DAEs1 0 CT -0 2陽+ %再+ 0.引再0 1 02珂X： -5%； -0 0 0%! + 2 4-1編寫Matlab代碼1. odefun=(t,x)-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);2. 2*x(1

33、)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);3. x(1)+x(2)+x(3)-1;% 微分方程組4. M=1 0 0;0 1 0;0 0 0;% 質(zhì)量矩陣5. options=odeset('mass',M);%對以 DAE 問題，mass屬性必須設(shè)置6. x0=0.8;0.1;0.1;% 初值7. t,x=ode15s(odefun,0 20,x0,options);% 這里好像不能使用 ode458. figure('numbertitle','off','name','DAE demo by Ma

34、tlabsky')plot(t,x)10.Iegend('x1(t)','x2(t)','x3(t)')復(fù)制代碼在介紹隱式ODEs時，我們說了 ode15i也可以求解DAEs原理就是將約束看成一個隱式 好，下面我們看看，驗證下所說的：1. odefun=(t,x,dx)dx(1)+0.2*x(1)-x(2)*x(3)-0.3*x(1)*x(2);2. dx(2)-2*x(1)*x(2)+5*x(2)*x(3)+2*x(2)A2;3. x(1)+x(2)+x(3)-1;4. %狀態(tài)變量初值，題目中給出5. x0=0.8;0.1;0.1;6.

35、 %該初值全部給定，按理說應(yīng)該全部為1，但是記住方程中有一個約束，故其實只有兩個獨立初值7. x0_fix=1;0;1;%隨意一個改為 0都可以，比如0;1;1或者1;1;08. %狀態(tài)變量一階微分初值，題目沒有給岀，可以任意寫dx0=1;1;1;9. %該初值一個都沒有給出，故全部為010. dxO_fix=O;O;O;11. %時間變量的初值12. tO=O;13. %計算相容初值14. xO,dxO=decic(odefun,O,xO,xO_fix,dxO,dxO_fix);15. xO,dxO % 相容初始條件16. res=ode15i(odefun,O,2O,xO,dxO);17.

36、 figure('numbertitle','off','name','DAE demo by Matlabsky')18. plot(res.x,res.y)其實，有些簡單DAEs是可以轉(zhuǎn)換為ODEs的，由于代數(shù)方程只是狀態(tài)變量的一個約束， 假如約束充分的話，我們就可以使用消參的思想，將約束反代回ODEs中，這樣就將那個約束方程消去了，最后只剩下那些微分方程了。對于本例，我們只要通過第三個方程表示出 再將x3反代回前兩個微分方程中則有' = -0.2 +Xj-兀2)+ 0. 3阿兀2花二 2 珂電5a(1 Aj a) 2

37、aj這個就是一個正常的ODEs了，好，下面我們就重新使用Matlab求解下看:1. x0=0.8;0.1;2. odefun=(t,x)-0.2*x(1)+x(2)*(1-x(1)-x(2)+0.3*x(1)*x(2)3. 2*x(1)*x(2)-5*x(2)*(1-x(1)-x(2)-2*x(2)*x(2);4. t,x=ode45(odefun,0,20,x0);5. figure('numbertitle','off','name','DAE demo by Matlabsky')6. plot(t,x)7. Iegend(

38、'x1(t)','x2(t)')教程延遲微分方程(DDE)的Matlab解法延遲微分方程的一般形式為：兀'(f) = 菇xQ-巧)/(一%),盂一耳)其中taui >0為狀態(tài)變量x(t)的延遲常數(shù)。Matlab中提供了求解這類方程的隱式Runge-Katta算法dde23()函數(shù)，可以直接求解延遲微分方程，該函數(shù)的格式為：sol = dde23(ddefun,tau,history,tspan,options )輸入?yún)?shù)：ddefun :描述延遲微分方程的句柄tau=tau1,taou, .taun:history :描述t < t0時的狀態(tài)

39、變量的值的函數(shù)，題目一般會直接給出tspan :同理 ode45注意返回參數(shù)sol.y是按行排列，只要與 ode*函數(shù)的轉(zhuǎn)置看下面的例子，進行說明。假設(shè)延遲微分方程組為+4®+2必)=4 欣) 其中& 莖 0時,= yl(£)= 0解(1)選取狀態(tài)變量。若想得出該方程的數(shù)值解，需要先將其變換成一階顯式微分方程組Aj(i)=A(/),可(d)二尹©), (/)=>'(/)(2) 計算每個狀態(tài)變量的一階微分珂'(Z) = 1 3x(f)疋2 ( 1 0.2 f (£ 0.5) Xj 0.5) 宅可O =碼)花'(f)二

40、4 嗎(£) 一 (f) 一 3 也(£)編寫延遲函數(shù)。從第一個方程可以看出x2延遲1, x1延遲0.5，我們?nèi)au=1 0.51. function dx=ddefun(t,x,z)2. %z(i,j)表示狀態(tài)變量xj延遲了tau(i)3. tau仁z(:,1);%第一列表示延遲了tau(1)=1的所有狀態(tài)變量4. tau2=z(:,2);5. %x2延遲了 1,故表示為z(2,1)或者tau1(2)6. %x1延遲了 0.5，故表示為z(1,2)或者tau2(1)7. dx=1-3*x(1)-tau1(2)-0.2*tau2(1F3-tau2(1)8. x(3)9.

41、4*x(1)-2*x(2)-3*x(3);(4)編寫主調(diào)函數(shù)1. tau=1 0.5;2. history=0 0 0'3. tspan=0 10;4. sol=dde23(ddefun,tau,history,tspan);5. figure('numbertitle','off','name','DDE demo by Matlabsky')6. plot(sol.x,sol.y(2,:)教程邊值問題(BVP)的 Matlab解法關(guān)于邊值微分方程問題相關(guān)連接有：【非線性微分方程邊值問題打靶算法】參見【線性微分方程邊值問

42、題打靶算法】參見【線性微分方程邊值問題有限差分算法】參見 【常微分方程 - 邊值問題六階精度解算器 (BVP6C)】 這里有好多 bvp4c 使用例子，不過 是英文的 到現(xiàn)在為止，前面講到的都是初值問題，也就是說 Matlab 的 ode* 系列解算器，默認將 tspan(1) 作為初值條件時的 t ，比如你將初值條件換為 x(2)=x '(2)=0 ， 那么 tspan(1) 就必 須是 2但是工程應(yīng)用中我們經(jīng)常遇到邊值問題， 這些是那些 ode* 函數(shù)無能為力的， 當(dāng)然我們可以 自己編寫函數(shù)求解 ( 比如 shooting) ，但是那個畢竟不是某些人能力所及的， 還好 Matlab

43、 中 提供了 bvp 解算器。solinit = bvpinit(x, yinit,params)sol = bvpsolver(odefun,bcfun,solinit,options )由于邊值問題可能有多解， 為了便于我們確定那個解是我們需要的， 所以必須使用 bvpinit 函數(shù)對初值進行估計解算器 (bvpsolver) ： Matlab 中提供了 bvp4c 和 bvp5c ，后者誤差控制更好些輸入?yún)?shù)：x ：需要計算的網(wǎng)格點，相當(dāng)于 ode* 的 tspanyinit ：狀態(tài)變量的初始猜測的值，可以是具體值，也可以是函數(shù)，類似與 ode* 的 x0params:微分方程中的其它未

44、知參數(shù)的初始猜測值odefun ：描述邊值問題微分方程的函數(shù)句柄bcfun :邊值函數(shù)，一般是雙邊值 (x 的上下限即認為兩個邊界 ) ，但也支持多邊值 (具體看幫 助)solinit :由 bvpinit 生成的初始化網(wǎng)格options : BVP解算器優(yōu)化參數(shù)，可以通過 bvpset設(shè)置，具體參數(shù)查看幫助輸出參數(shù)：大部分同理ode45，這里不詳細介紹了入也是未好下面我們以例子說事，考察包含未知參數(shù)入的二階微分方程，求解入，由于知的(這就是在params中說到的其他位置參數(shù))，故需要 三個邊值條件y"+ (2 - 2q cos 2x) = 0 其中q = 5,y0) =y'

45、(7f) = 0, y (0) = 1解(1)選用狀態(tài)變量(2)求每個狀態(tài)變量的一階微分' = y'= x3<花 1 =丸 + 2q cos 2： j(3) 編程實現(xiàn)求解1. q=5;2. lambda = 15;%未知參數(shù)猜測值3. x=linspace(0,pi,10);% 需要計算的點4. yinit=(x)cos(4*x);-4*sin(4*x);% 使用函數(shù)估計目標(biāo)值5. solinit = bvpinit(x,yinit,lambda);% 生成計算網(wǎng)格6. odefun=(x,y,lambda)y(2);-(lambda - 2*q*cos(2*x)*y(1

46、);% 邊值微分方程7. bcfun=(ya,yb,lamdba)ya(2);yb(2);ya(1)-1;% 邊界條件8. sol = bvp4c(odefun,bcfun,solinit);9. fprintf('The fourth eigenvalue is approximately %7.3f.n',sol.parameters)10. xint = linspace(0,pi);11. Sxint = deval(sol,xint);% 計算 xint 對應(yīng)的 y 值12. figure('name','BVP Demo by Matlabs

47、ky', ' numbertitle ' , ' off ')13. plot(xint,Sxint(1,:)% 繪制 xint 和 x1=y 的圖形，即 x 與 y14. axis(0 pi -1 1.1)15. title('Eigenfunction of Mathieu"s equation.')16. xlabel('x')17. ylabel('solution y')已;un lion of弓“川前沖般偏微分方程組(PDEs)的MATLAB求解一、pdepe()函數(shù)說明MATLAB語

48、言提供了 pdepe()函數(shù)，可以直接求解一般偏微分方程(組)，它的調(diào)用格式為1. sol=pdepe(m,pdefun,pdeic,pdebc,x,t)【輸入?yún)?shù)】pdefun：是PDE的問題描述函數(shù)，它必須換成下面的標(biāo)準形式0 （和浮）申=計（兀覚見舟）十衛(wèi),申）（式1）辦 atoxaxdx這樣，PDE就可以編寫下面的入口函數(shù)1. c,f,s=pdefun（x,t,u,du）c,f,s這三個函m,x,t就是對應(yīng)于（式1）中相關(guān)參數(shù)，du是u的一階導(dǎo)數(shù)，由給定的輸入變量即可表示出出 數(shù)pdebc：是PDE的邊界條件描述函數(shù)，必須先化為下面的形式戸（心3 +0（和忌）護r） = 0于是邊值條件

49、可以編寫下面函數(shù)描述為1. pa,qa,pb,qb=pdebc（x,t,u,du）其中a表示下邊界，b表示下邊界pdeic：是PDE的初值條件，必須化為下面的形式股我們使用下面的簡單的函數(shù)來描述為1. uO=pdeic(x)m,x,t:就是對應(yīng)于(式1)中相關(guān)參數(shù)【輸出參數(shù)】sol：是一個三維數(shù)組，sol(:,:,i)表示ui的解，換句話說 uk對應(yīng)x(i)和t(j)時的解為sol(i,j,k)通過sol，我們可以使用pdeval()直接計算某個點的函數(shù)值1、PDEs問題實例講解試求解下面的偏微分魯"024勢-叫-嗎) 竽"17券+盹宀)其中,F(x)=嚴-11.46r且滿足初始條件 ':-'"，及邊界條件學(xué)3丿)=眄(6Q = o.坷(1丿)=h學(xué)1丿)=oOXOX【解】（1 ）對照給出的偏微分方程，根據(jù)標(biāo)注形式，則原方程可以改寫為（旳一旳） F （叭一妁）注意：1;1那里是“.*而不是“，”表示對應(yīng)元素相乘，而不是矩陣的相乘，不要問為什