MATLAB 程式設計/高階主題/數值運算/常微分方程

所有微分方程都具有相同的語法，您必須使用該語法，以及相同的輸入和輸出引數。所有求解器都需要三個輸入引數：您要解決的微分方程的函式控制代碼、要積分的時間間隔和一組初始條件。假設我們要解決簡單的微分方程 y' = y，y(0) = 1，其真解為 y(t) = e^t。假設我們希望在 t = 0 到 t = 1 之間得到解。

要使用函式控制代碼，您必須首先建立一個包含該函式的 M 檔案，如下所示

function Yprime = func(t, y)
Yprime = 0;
Yprime = y;

請注意，即使時間引數未在微分方程中使用，您也必須包含它。將 Yprime 初始化為一個零陣列將在您嘗試解決多個函式時為您節省麻煩；Yprime 必須作為垂直陣列返回，但是，如果您沒有將其初始化為垂直陣列（或在最後進行轉置），它預設情況下會返回一個水平陣列。

還要注意，除了 ode15i 之外，必須顯式地針對 y' 求解該函式。

建立此檔案後，按順序（函式、時間間隔、初始條件）使用引數呼叫 ODE 函式

>> ode45(@func, [0,1], 1)

另一種方法是使用匿名函式，這僅在您只有一個函式時有用（否則您必須使用函式控制代碼）。您還必須根據因變數和時間來宣告匿名函式

>> func = @(t,y) y;
>> ode45(func, [0,1], 1)

在沒有輸入引數的情況下呼叫 ODE 函式會生成解的圖形。使用一個輸出引數呼叫它會返回一個結構陣列

>> Struct = ode45(func, [0,1], 1)
 Struct = solver: 'ode45'
          extdata: [1x1 struct]
          x: [1x11 double]
          y: [1x11 double]
          stats: [1x1 struct]
          idata: [1x1 struct]

此資料主要用於將來呼叫“deval”函式以在您想要的時間獲取答案

>> Solution = deval(Struct, 0.5)
 Solution = 1.6487

Deval 還可以透過包含第二個輸出引數來返回感興趣點的導數。

>> [Solution, Derivative] = deval(Struct, 0.5)
Solution = 1.6487
Derivative = 1.6487

由於 e^x 的導數是其本身，因此導數和解在這裡相同是有道理的。

使用兩個輸出引數呼叫 ode45 會返回兩個資料列表；首先是 t，然後是適當大小矩陣中的自變數。

ODE 選項

MATLAB 為您提供了相當多的選項來修改它求解微分方程的方式。幫助檔案對所有選項都做了很好的描述，因此這裡將不再詳細介紹。要獲取列表，請使用幫助函式

>>help odeset

要獲取不同選項名稱的列表以及要傳遞給它的內容，只需在命令提示符中鍵入“odeset”。它返回一個數據型別或一個有限的選項列表。它還以括號形式列出了所有選項的預設值。

要設定特定選項或選項列表，請先鍵入選項名稱，然後鍵入您想要的選項值。例如，假設您要將預設相對容差從 10^-3 調整為 10^-4。您將按如下方式呼叫“odeset”

>> option = odeset('RelTol', 10^-4);

請注意，選項名稱必須作為字串傳遞，否則您很可能會收到“未定義變數”錯誤。可以透過將所有選項都放在一個列表中來一次傳遞多個選項

>> option = odeset('RelTol', 10^-4, 'AbsTol', 10^-7, 'MassSingular', 'no');

然後可以將選項結構作為第四個（可選）輸入引數傳遞給 ode 函式

>> [T,X] = ode45(@func, [0,1], 1, option);

這將返回比預設值更準確的值，因為誤差容差更嚴格。它還應更快地計算，因為 MATLAB 不會檢查這是否是微分代數方程（這就是 MassSinglular 選項的作用；它通常設定為“maybe”，因此 MATLAB 會自行檢查）。

ODE 求解器

MATLAB 並不只有一種 ODE 求解器，截至 MATLAB 7.1 版本，它有八種。有些求解器更適合解決某些問題，而不是其他問題，這就是為什麼所有這些求解器都被包含在內的原因。MATLAB 幫助檔案中列出了每個函式可以執行的操作，但這裡提供一個快速摘要，大致按照您應該嘗試它們的順序（除非您已經知道問題的分類）

大多數問題都可以使用 ode45 求解，由於這是速度和精度的最佳折衷方案，因此它應該是您首先使用的求解器。

如果您需要非常嚴格的誤差容差或大量資料點，請使用 ode113。

如果您需要相對寬鬆的誤差容差或問題在 ode45 中很慢，請嘗試使用 ode23。

如果問題確實是剛性的，而 ode45 失敗，請對於寬鬆的容差使用 ode23tb，對於稍微嚴格的容差（具有恆定質量矩陣）使用 ode23s，對於更嚴格的容差或非恆定質量矩陣使用 ode15s。

如果您需要在剛性問題上進行數值阻尼的解，請使用 ode23t。

s 表示該演算法適用於剛性問題。

還有另一種 ODE 求解器，它是特殊的

只能使用 ode15i 求解隱式問題。

由於 ode15i 在語法方面有一些區別，因此將在下一節中進行討論。

隱式 ODE：ODE15i

由於 ode15i 是唯一一個求解隱式方程的 ODE 求解器，因此它在如何輸入函式方面必須有一些特殊的語法規則。

如果您使用的是 ode15i，請按如下方式宣告該函式

function Z = func(t,y,Yprime)
Z = 0;
Z = y - Yprime;

請注意，對於 ode15i，您必須將函式置於標準形式（針對 0 求解），而對於所有其他 ODE 函式，您必須顯式地針對 y' 求解。還要注意，您必須宣告三個輸入引數，而不是通常的兩個引數。

當您呼叫 ode15i 時，您不僅必須包含 y 的初始條件，還必須包含 Yprime 的初始條件。Yprime 的初始條件進入第四個引數

>> [t,x] = ode15i(@func, [0,1], 1, 1);

這將返回與 ode45 相似的結果（用於顯式方程），但資料點更少。

選項結構作為可選的第五個引數傳遞給 ode15i。ode15i 的所有輸出選項與其他 ODE 求解器相同。