在工程實踐中,我們經常遇到一些odes,其中某些解變換緩慢,另一些變化很快,且相差懸殊的微分方程,這就是所謂的剛性問題(stiff),對於所有解的變化相當我們則稱為非剛性問題(nonstiff)。
變步長模式解法器有:ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb和discrete。
a) ode45:預設值,四/五階龍格-庫塔法,適用於大多數連續或離散系統,但不適用於剛性(stiff)系統。它是單步解法器,也就是,在計算y(tn)時,它僅需要最近處理時刻的結果y(tn-1)。一般來說,面對乙個**問題最好是首先試試ode45。
b) ode23:二/三階龍格-庫塔法,它在誤差限要求不高和求解的問題不太難的情況下,可能會比ode45更有效。也是乙個單步解法器。
c) ode113:是一種階數可變的解法器,它在誤差容許要求嚴格的情況下通常比ode45有效。ode113是一種多步解法器,也就是在計算當前時刻輸出時,它需要以前多個時刻的解。
d) ode15s:是一種基於數字微分公式的解法器(ndfs)。也是一種多步解法器。適用於剛性系統,當使用者估計要解決的問題是比較困難的,或者不能使用ode45,或者即使使用效果也不好,就可以用ode15s。
e) ode23s:它是一種單步解法器,專門應用於剛性系統,在弱誤差允許下的效果好於ode15s。它能解決某些ode15s所不能有效解決的stiff問題。
f) ode23t:是梯形規則的一種自由插值實現。這種解法器適用於求解適度stiff的問題而使用者又需要乙個無數字振盪的解法器的情況。
g)ode23tb:是tr-bdf2的一種實現, tr-bdf2 是具有兩個階段的隱式龍格-庫塔公式。
h)discrtet:當simulink檢查到模型沒有連續狀態時使用它。
固定步長模式解法器有:ode5,ode4,ode3,ode2,ode1和discrete。
a) ode5:預設值,是ode45的固定步長版本,適用於大多數連續或離散系統,不適用於剛性系統。
b) ode4:四階龍格-庫塔法,具有一定的計算精度。
c) ode3:固定步長的二/三階龍格-庫塔法。
d) ode2:改進的尤拉法。
e) ode1:尤拉法。
f) discrete:是乙個實現積分的固定步長解法器,它適合於離散無連續狀態的系統。
表2-3
函式指令
含 義
函 數
含 義
求解器solver
ode23
普通2-3階法解ode
odefile
包含ode的檔案
ode23s
低階法解剛性ode
選項odeset
建立、更改solver選項
ode23t
解適度剛性ode
odeget
讀取solver的設定值
ode23tb
低階法解剛性ode
輸出odeplot
ode的時間序列圖
ode45
普通4-5階法解ode
odephas2
ode的二維相平面圖
ode15s
變階法解剛性ode
odephas3
ode的三維相平面圖
ode113
普通變階法解ode
odeprint
在命令視窗輸出結果
3.因為沒有一種演算法可以有效地解決所有的ode問題,為此,matlab提供了多種求解器solver,對於不同的ode問題,採用不同的solver。
表2-4 不同求解器solver的特點
求解器solver
ode型別
特點說明
ode45
非剛性一步演算法;4,5階runge-kutta方程;累計截斷誤差達(△x)3
大部分場合的首選演算法
ode23
非剛性一步演算法;2,3階runge-kutta方程;累計截斷誤差達(△x)3
使用於精度較低的情形
ode113
非剛性多步法;adams演算法;高低精度均可到10-3~10-6
計算時間比ode45短
ode23t
適度剛性
採用梯形演算法
適度剛性情形
ode15s
剛性多步法;gear』s反向數值微分;精度中等
若ode45失效時,可嘗試使用
ode23s
剛性一步法;2階rosebrock演算法;低精度
當精度較低時,計算時間比ode15s短
ode23tb
剛性梯形演算法;低精度
當精度較低時,計算時間比ode15s短
4.在計算過程中,使用者可以對求解指令solver中的具體執行引數進行設定(如絕對誤差、相對誤差、步長等)。
表2-5 solver中options的屬性
屬性名取值
含義abstol
有效值:正實數或向量
預設值:1e-6
絕對誤差對應於解向量中的所有元素;向量則分別對應於解向量中的每一分量
reltol
有效值:正實數
預設值:1e-3
相對誤差對應於解向量中的所有元素。在每步(第k步)計算過程中,誤差估計為:
e(k)<=max(reltol*abs(y(k)),abstol(k))
normcontrol
有效值:on、off
預設值:off
為『on』時,控制解向量範數的相對誤差,使每步計算中,滿足:norm(e)<=max(reltol*norm(y),abstol)
events
有效值:on、off
為『on』時,返回相應的事件記錄
outputfcn
有效值:odeplot、odephas2、odephas3、odeprint
預設值:odeplot
若無輸出參量,則solver將執行下面操作之一:
畫出解向量中各元素隨時間的變化;
畫出解向量中前兩個分量構成的相平面圖;
畫出解向量中前三個分量構成的三維相空間圖;
隨計算過程,顯示解向量
outputsel
有效值:正整數向量
預設值:
若不使用預設設定,則outputfcn所表現的是那些正整數指定的解向量中的分量的曲線或資料。若為預設值時,則預設地按上面情形進行操作
refine
有效值:正整數k>1
預設值:k = 1
若k>1,則增加每個積分步中的資料點記錄,使解曲線更加的光滑
jacobian
有效值:on、off
預設值:off
若為『on』時,返回相應的ode函式的jacobi矩陣
jpattern
有效值:on、off
預設值:off
為『on』時,返回相應的ode函式的稀疏jacobi矩陣
mass
有效值:none、m、
m(t)、m(t,y)
預設值:none
m:不隨時間變化的常數矩陣
m(t):隨時間變化的矩陣
m(t,y):隨時間、地點變化的矩陣
maxstep
有效值:正實數
預設值:tspans/10
最大積分步長
y(0)=1,y』(0)=0
令x1=y,x2=dy/dx,則
dx1/dt = x2
dx2/dt = μ(1-x2)-x1
編寫函式檔案verderpol.m:
function xprime = verderpol(t,x)
global mu
xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)];
再在命令視窗中執行:
>>global mu
>>mu = 7;
>>y0=[1;0]
>>[t,x] = ode45(『verderpol』,0,40,y0);
>>x1=x(:,1);x2=x(:,2);
>>plot(t,x1,t,x2)
圖形結果為圖2-20。
用Matlab實現簡單有限體積求解器
用於瞬態對流擴散pde的簡單而通用的fvm求解器 乙個簡單的有限體積工具 這段 是化學 石油工程師開發一種簡單的工具來求解對流擴散方程的一般形式的結果 general equationt u d 在1d,1d軸對稱 徑向 2d,2d軸對稱 圓柱 和3 d域上的簡單均勻 不均勻網格上。該 在整個或部分...
matlab求解振動方程
看了一篇柱塞幫浦離散化動力學建模的文章,感覺還挺有意思,於是嘗試做一下 二 matlab下的動力學方程總結 斜盤式軸向柱塞幫浦是一類常見的柱塞幫浦,本文以 型斜盤式軸向柱塞幫浦為研究物件,研究幫浦內機械振動的傳遞問題。由於該幫浦傳動軸與缸體之間為過盈配合,且柱塞滑靴元件位於缸體的柱塞腔內,因此,將傳...
MATLAB求解矩陣函式
3.一般矩陣運算函式不可用 4.矩陣函式求解函式funm 5.多多點贊關注,多多交流 這裡的 1 2 節涉及到矩陣函式的一般求法。第 3 節演示了為什麼不能用普通的運算函式求矩陣函式。第 4 節介紹了matlab內建的矩陣函式求解函式。如果趕時間,直接看第 4 節即可!總的來說,矩陣函式的求解方式和...