4.5　數值計算

數值計算是指利用計算機求數學問題（例如，函數的零點、極值、積分和微分以及微分方程）近似解的方法。常用的數值分析有求函數的最小值、求過零點、數值微分、數值積分和解微分方程等。

4.5.1　函數極值

數學上利用計算函數的導數來確定函數的最大值點和最小值點，然而，很多函數很難找到導數為零的點。為此，可以通過數值分析來確定函數的極值點。MATLAB只有處理極小值的函數，沒有專門求極大值的函數，因為f（x）的極大值問題等價于﹣f（x）的極小值問題。MATLAB求函數的極小值使用fminbnd和fminsearch函數。

1．一元函數的極值

fminbnd函數可以獲得一元函數在給定區間內的最小值，函數調用格式如下：

其中，fun是函數的句柄或匿名函數；x1和x2是尋找函數最小值的區間范圍（x1＜x＜x2）；x為在給定區間內，極值所在的橫坐標。

其中，y為求得的函數極值點處的函數值。

【例4-16】　已知y=e^﹣0.2xsin（x），在0≤x≤5π區間內，使用fminbnd函數獲取y函數的極小值。

程序代碼如下：

程序運行結果如下，圖4-7是函數在區間[0，5?pi]的函數曲線圖。

由圖4-7可知，函數在x=4.5附近出現極小值點，極小值約為﹣0.4，驗證了用極小值fminbnd函數求的極小值點和極小值是正確的。

2．多元函數的極值

fminsearch函數可以獲得多元函數的最小值，使用該函數時需要指定開始的初始值，獲得初始值附近的局部最小值。該函數的調用格式如下：

其中，fun是多元函數的句柄或匿名函數；x0是給定的初始值；x是最小值的取值點；y是返回的最小值，可以省略。

【例4-17】　使用fminsearch函數獲取f（x，y）二元函數在初始值（0，0）附近的極小值，已知f（x，y）=100（y﹣x²）²＋（1﹣x）²。

程序代碼如下：

程序運行結果如下：

由結果可知，由函數fminsearch計算出局部最小值點是[1，1]，最小值為y1=3.6862e﹣10，和理論上是一致的。

4.5.2　函數零點

一元函數f（x）的過零點的求解相當于求解f（x）=0方程的根，MATLAB可以使用fzero函數實現，需要指定一個初始值，在初始值附近查找函數值變號時的過零點，也可以根據指定區間來求過零點。該函數的調用格式為

其中，x為過零點的位置，如果找不到，則返回NaN；y是指函數在零點處函數的值；fun是函數句柄或匿名函數；x0是一個初始值或初始值區間。

需要指出，fzero函數只能返回一個局部零點，不能找出所有的零點，因此需要設定零點的范圍。

【例4-18】　使用fzero函數求f（x）=x²﹣5x＋4分別在初始值x=0和x=5附近的過零點，并求出過零點函數的值。

程序代碼如下：

程序運行結果如下：

由結果可知，用fzero函數可以求在初始值x0附近的函數過零點。不同的零點，需要設置不同的初始值x0。

4.5.3　數值差分

任意函數f（x）在x點的前向差分定義為

稱Δf（x）為函數f（x）在x點處以h（h＞0）為步長的前向差分。

在MATLAB中，沒有直接求數值導數的函數，只有計算前向差分的函數diff，其調用格式為

例如，已知矩陣，分別求矩陣A行和列的一階和二階前向差分。

4.5.4　數值積分

數值積分是研究定積分的數值求解的方法。MATLAB提供了很多種求數值積分的函數，主要包括一重積分和二重積分兩類函數。

1．一重數值積分

MATLAB提供了quad函數和quadl函數求一重積分。它們的調用格式分別如下：

它是一種采用自適應的Simpson方法的一重數值積分，其中，fun為被積函數，函數句柄；a和b為定積分的下限和上限；tol為絕對誤差容限值，默認是10^﹣6；trace控制是否展現積分過程，當trace取非0，則展現積分過程，默認取0。

它是一種采用自適應的Lobatto方法的一重數值積分，參數定義和quad一樣。

【例4-19】　分別使用quad函數和quadl函數求的數值積分。

程序代碼如下：

程序運行結果如下：

其中，迭代過程最后一列的和為數值積分q3的值。

2．多重數值積分

MATLAB提供了dblquad函數和triplequad函數求二重積分和三重積分。它們的調用格式如下：

函數的參數定義和一重積分一樣。

例如，求二重數值積分。

代碼如下：

4.5.5　常微分方程求解

MATLAB為解常微分方程提供了多種數值求解方法，包括ode45、ode23、ode113、ode15s、ode23s、ode23t和ode23tb等函數，用得最多的是4/5階龍格－庫塔法ode45函數。該函數的調用格式如下：

其中，fun是待解微分方程的函數句柄；ts是自變量范圍，可以是范圍[t0，tf]，也可以是向量[t0，…，tf]；y0是初始值，y0和y具有相同長度的列向量；options是設定微分方程解法器的參數，可以省略，也可以由odeset函數獲得。

需要注意，用ode45求解時，需要將高階微分方程y^（n）=f（t，y，y′，…，y^（n﹣1）），改寫為一階微分方程組，通常解法是，假設y₁=y，從而y₁=y，y₂=y′，…，y_n=y^（n﹣1），于是高階微分方程可以轉換為下述常微分方程組求解：

【例4-20】　已知二階微分方程，試用ode45函數解微分方程，作出y—t的關系曲線圖。

（1）首先把二階微分方程改寫為一階微分方程組。

令y₁=y，，則

（2）程序代碼如下：

程序運行結果如圖4-8所示。