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所示。