ODE45函數(shù)的使用——翻譯.doc

在Matlab中使用ode45簡介Matlab中常微分方程常用的函數(shù)是ODE45，這個(gè)函數(shù)能夠利用-龍哥庫塔法-有效求解帶時(shí)間變量步長的計(jì)算。Ode45用于求解如下的一般問題：（1）其中，時(shí)間t是獨(dú)立變量，x為時(shí)間相關(guān)矢量，是時(shí)間t和x的函數(shù)。當(dāng)（1）右邊的是固定的，且給定x的初始值，那么問題的解是唯一的。在ME175中，解法是不完整的，但是只要你解決了問題，就可以獲得ODE代表的系統(tǒng)運(yùn)動(dòng)趨勢(shì)。這有利于得到一個(gè)直觀的印象，看起來很復(fù)雜的常微分方程，代表的質(zhì)點(diǎn)運(yùn)動(dòng)軌跡確實(shí)簡單明了的。以下簡要解釋如何得到運(yùn)動(dòng)軌跡：第一步：對(duì)給定的ODE方程進(jìn)行降階處理，得到一系列一階方程這就是你要做的第一步，在一張草稿紙上處理。例如，給定ODE方程如下：（2）對(duì)本問題，矢量x有兩個(gè)組成分量：y和，或（3）且（4）其中，用（3）中的式子代表了y，于是把（2）改寫為（4）。如果求解的問題有更多階數(shù)更多變量呢？例如，我們除了有上面的方程（2），同時(shí)還有以下的方程：（5）那么，我們可以通過構(gòu)造更大的矢量x同時(shí)求解y，z： （6）然后 （7）以及（8）其中，y變量和z變量的放置位置對(duì)求解不造成影響。實(shí)際上，任意次序都是有效的，例如 和但是重要的是，在整個(gè)計(jì)算過程中，你使用的順序都必須和一階ODE方程中定義的變量順序相同。之后，如果你使用的是（7）中給定的的式子，那么系統(tǒng)的一階ODE方程，由以下方程組組成。（10）而涉及的表征變量結(jié)果如下：（11）基本上，可以處理任意數(shù)量的高階ODE方程。重要的是把它們處理成多個(gè)一階的ODE方程，并且確保記住被求解的矢量X中，不同變量所分配的順序。第二步 編寫代碼既然你已經(jīng)有所求解問題的一階格式，在你編程的主要代碼中，將會(huì)用到以下的命令fname是函數(shù)的M文件名用于求解方程（1）右邊代數(shù)式的值。這個(gè)函數(shù)將被輸入一階ODE系統(tǒng)中，并且被積分（見（10），（11）。后面，將會(huì)更詳細(xì)的解釋。注：當(dāng)然關(guān)于ODE45如何積分給定的方程有細(xì)微的差別，但是對(duì)于簡單的問題，不分先后次序的積分，是可以接受的。tspan 是矢量定義了積分的起始點(diǎn)和終點(diǎn)，同時(shí)也定義了時(shí)間步長。例如，我們需要積分t=0到t=10，希望步數(shù)是100步，那么tspan=0:0.1:10或者tspan=linspace（0,10,100).xinit是初始條件矢量。確保初始值的順序和給定的x中變量和它倒數(shù)的順序是一致的。同時(shí)注意如果x有5個(gè)變量，那么同時(shí)要輸入5個(gè)初始值。option這個(gè)在matlab的幫助文件中有很好的說明。對(duì)于大部分的問題，使用默認(rèn)值就可以滿足計(jì)算要求。t是獨(dú)立變量，計(jì)算數(shù)組x在時(shí)間點(diǎn)t的數(shù)值。這個(gè)矢量不必等于tspan，ODE45自動(dòng)調(diào)節(jié)步數(shù)以取得最大的效率和精確度。（在快速變化部分采用小步長，在變化緩慢部分采用大步長）。x相關(guān)內(nèi)容如下。X是數(shù)組或矩陣，大小為length（t）*length（xinit）。每一列x代表不同的因變量。例如，為簡單假定t=0,1,2.,10,將會(huì)計(jì)算函數(shù)在11個(gè)點(diǎn)的值。（12）如果是x的第四個(gè)變量，那么得到了在t=0時(shí)候的值，得到了在t=6時(shí)候的值，得到了t=10時(shí)候的值。簡而言之，代表x的第k個(gè)變量，k=1與變量y相關(guān)，k=2與變量相關(guān)。計(jì)算所有變量在某一時(shí)間點(diǎn)j的數(shù)值注：在產(chǎn)生hokey pokey舞蹈前，史前兒童圍坐在篝火前齊唱：You put your left foot inYou put your left foot outYou put your left foot in And you shake it all about當(dāng)你使用matlab函數(shù)ODE45及時(shí)完成作業(yè)時(shí)，x就是要做得全部內(nèi)容。不幸的是由于缺乏matlab軟件，使得這本書過時(shí)了。命令的作用是重新定義變量。、，然后，如果使用變量順序?yàn)椋瑧?yīng)該這么寫程序：當(dāng)然，也不應(yīng)該認(rèn)為定義y，ydot麻煩。直接表達(dá)為x的形式（例如，使用代表y），清晰的定義方式有利于后面的調(diào)試。以下，你將以（或是被要求）繪圖的形式描述感興趣的軌跡：質(zhì)點(diǎn)隨時(shí)間運(yùn)動(dòng)軌跡，在平面中表示角度和徑向關(guān)系等。繪圖和繪制子圖的命令在matlab幫助文件中有清晰的說明，這里不再詳細(xì)說明。記住如果你想在一張圖中放多個(gè)圖，應(yīng)該使用子圖的概念，當(dāng)然在一張圖片中畫多個(gè)圖，不是一個(gè)好主意。別忘記給圖加標(biāo)簽：包括標(biāo)題，x軸，y軸的含義，如果多條曲線應(yīng)該分別標(biāo)明。最后，請(qǐng)注意在中包含的僅僅是變量而已，依據(jù)自己的喜好使用字母，T替換t，x0替換xinit都是可以的。只有記住使用新變量名，之后的每個(gè)引用都用一樣的名稱。另一個(gè)普遍的錯(cuò)誤在于，同一變量的重復(fù)定義。例如，定義，如果足夠幸運(yùn)的話，會(huì)有錯(cuò)誤警告；不幸的話，這種錯(cuò)誤很難發(fā)現(xiàn)，要花數(shù)小時(shí)時(shí)間檢查您的程序以解決問題。此外，fname是什么呢？回憶下，我們還沒有告訴matlab程序應(yīng)該對(duì)什么函數(shù)進(jìn)行積分，是吧？這就是為什么需要fname文件，fname文件含有所有之前在稿紙上重寫的ODE一階函數(shù)。你可以對(duì)這個(gè)文件起任意的名稱，只要與中使用的fname一致。例如，你對(duì)fname取名superman那么是對(duì)的，而就不正確了。更進(jìn)一步說，函數(shù)不必須像在原始代碼中寫的在同一個(gè)m文件中，一些人喜歡在程序末尾書寫子程序，特別是代碼不長，比較簡單的時(shí)候。例如，你的代碼名稱ME175example文件，那么m文件將如下：Function dxdt=ME175example(t,x)%這里 t,x 和 dxdt are 只是變量而已。你可以起任意名稱%只要 t 是獨(dú)立變量，而X是因變量%dxdt 是推到的一階因變量% 定義常數(shù)m= 1;%定義變量使之清晰易懂%Recall that x = y, ydot, z, zdot, zdotdoty=x(1);ydot=x(2);z=x(3);zdot=x(4);zdotdot=x(5);%注意x僅僅是1列五行的數(shù)組%t,x = ode45(fname, tspan, xinit, options)%數(shù)組dxdt與x的大小相同dxdt = zeros(size(x);dxdt(1) = ydot;dxdt(2) = 1/m(5 x(2)exp(y) +y2); %This is ydotdotdxdt(3) = zdot;dxdt(4) = zdotdot;dxdt(5) = t-zdotdot+sin(z); %This is zdotdotdot%Note that the input arguments must be t and x (in that order) even in the case where t is notexplicitly used in the function.基本模板以下是基本模板，當(dāng)你想對(duì)一個(gè)高階常微分方程進(jìn)行積分時(shí)，把它復(fù)制黏貼到Matlab中Function 任何你想要的名字%定義起始時(shí)間tstat，終止時(shí)間tend，時(shí)間步數(shù)nTstart=？;Tend=？;N=?;Tspan=linspace(Tstart, Tend, N);%定義初始值，確保正確的順序Xinit=.; .; .; .;.;%獲得矢量x。把option設(shè)置為默認(rèn)值即可%定義輸出變量.%所需要畫圖的函數(shù) Subplot（?, ?, ?)Plot(需要的圖像）%plot3（）畫3D圖Title（ )Xlabel( )Ylabel( )Zlabel( )-%積分函數(shù)，可以作為獨(dú)立m文件，或在本程序底部，語法如下Function dxdt=積分函數(shù)（t，x）%定義積分中使用的常數(shù)%定義新變量；%也可以不定義，這樣的好處是，便于其他人閱讀你的程序%寫下你得到的一階ODE方程dxdt=zeros（size（x）dxdt（1）=?;dxdt（2）=?;dxdt（3）=?;Etc.小結(jié)：如果所有的步驟都對(duì)了，方程也是正確的。你將發(fā)現(xiàn)你的到的關(guān)系圖是很優(yōu)美的，不用你艱難的推到。干的好！最后，同樣地問題又多種不同的解法。多用matlab嘗試，你將有所發(fā)現(xiàn)，比起其他算法有些算法非常有效率。我使用Matlab三年了，