第13章数值分析 每当难以对一个函数进行积分、微分或者解析上确定一些特殊的值时,就可以借助计 算机在数值上近似所需的结果。这在计算机科学和数学领域,称之为数值分析。至此,可以 猜到, MATLAB提供了解决这些问题的工具。本章将介绍这些工具的使用 13.1绘图 说到绘图,只要计算函数在某一区间的值,并且画出结果向量,这样就得到了函数的 图形。在大多数情况下,这就足够了。然而,有时一个函数在某一区间是平坦的并且无激励 而在其它区间却失控。在这种情况下,运用传统的绘图方法会导致图形与函数真正的特性相 去甚远。 MATLAB提供了一个称为 fplot的巧妙的绘图函数。该函数细致地计算要绘图的函 数,并且确保在输出的图形中表示出所有的奇异点。该函数的输入需要知道以字符串表示的 被画函数的名称以及2元素数组表示的绘图区间。例如: >>fplot(humps,0 2) >>title(FPLOT OF HUMPS) 在0和2之间计算函数 humps,并显示该函数的图形。(见图13.1)。 FPLOT OF HUMPS 100 0.5 图13.1函数 humps的图形
第 13 章 数值分析 每当难以对一个函数进行积分、微分或者解析上确定一些特殊的值时，就可以借助计 算机在数值上近似所需的结果。这在计算机科学和数学领域，称之为数值分析。至此，可以 猜到，MATLAB 提供了解决这些问题的工具。本章将介绍这些工具的使用。 13.1 绘图 说到绘图，只要计算函数在某一区间的值，并且画出结果向量，这样就得到了函数的 图形。在大多数情况下，这就足够了。然而，有时一个函数在某一区间是平坦的并且无激励， 而在其它区间却失控。在这种情况下，运用传统的绘图方法会导致图形与函数真正的特性相 去甚远。MATLAB 提供了一个称为 fplot 的巧妙的绘图函数。该函数细致地计算要绘图的函 数，并且确保在输出的图形中表示出所有的奇异点。该函数的输入需要知道以字符串表示的 被画函数的名称以及 2 元素数组表示的绘图区间。例如： >>fplot(‘ humps ‘ , [0 2]) >>title(‘ FPLOT OF HUMPS ‘) 在 0 和 2 之间计算函数 humps，并显示该函数的图形。（见图 13.1)。 0 0.5 1 1.5 2 -20 0 20 40 60 80 100 FPLOT OF HUMPS 图 13.1 函数 humps 的图形
在这个例子中, humps·是 MATLAB的M文件函数 function y=humps(x) HUMPS A function used by QUADDEMO, ZERODEMO and FPLOTDEMO o HUMPS(X) is a function with strong maxima near x= 3 and x.9 See QUADDEMO, ZERODEMO and FPLOTDEMO Copyright(c)1984-93 by The Math Works, Inc y=1/(x-.3).^2+.01)+1,/(x-.9)^2+.04)-6; fplot适用于任何具有单输入和单输出向量的函数M文件。即如同 humps,输出变量y 返回一个与输入x同样大小的数组,在数组到数组意义上y和x有联系。在使用 fplot(以 及其它数值分析函数)的过程中,最普遍犯的错误是忘记把函数名加上引号。即 fplot需要 知道字符串形式的函数名。如果输入 fplot( humps,|0,2), MATLAB认为hmps是工作空 间中的一个变量,而不是函数的名称。注意把变量 humps定义为所需要的字符串,就可避 免这个问题。 >>humps="humps >>fplot(hump, [0 这时, MATLAB从变量 humps中获得字符串 humps 对于可表示成一个字符串的简单的函数,如y=2e-sin(x), fplot绘制这类函数的曲 线时,不用建立M文件,只需把x当作自变量,把被绘图的函数写成一个完整的字符串 >>f=2*exp(-x).*sin(x) 式中,运用数组乘法定义了函数f(x)=2e-sin(x) >> fplot(f,[08]) >>title(f), xlabel(x)
在这个例子中，‘ humps ‘是 MATLAB 的 M 文件函数。 function y=humps(x) % HUMPS A function used by QUADDEMO, ZERODEMO and FPLOTDEMO. % HUMPS(X) is a function with strong maxima near x= .3 and x= .9. % See QUADDEMO, ZERODEMO and FPLOTDEMO. % Copyright (c) 1984-93 by The MathWorks, Inc. y=1 ./ ((x - .3) .^ 2+ .01)+1 ./ ((x - .9) .^ 2+ .04) - 6; fplot 适用于任何具有单输入和单输出向量的函数 M 文件。即如同 humps，输出变量 y 返回一个与输入 x 同样大小的数组，在数组到数组意义上 y 和 x 有联系。在使用 fplot（以 及其它数值分析函数）的过程中，最普遍犯的错误是忘记把函数名加上引号。即 fplot 需要 知道字符串形式的函数名。如果输入 fplot(humps , [0 , 2])，MATLAB 认为 humps 是工作空 间中的一个变量，而不是函数的名称。注意把变量 humps 定义为所需要的字符串，就可避 免这个问题。 >>humps=‘ humps ‘; >>fplot(hump , [0 2]) 这时，MATLAB 从变量 humps 中获得字符串‘ humps‘。 对于可表示成一个字符串的简单的函数，如 y e x x = − 2 sin( ) ，fplot 绘制这类函数的曲 线时，不用建立 M 文件，只需把 x 当作自变量，把被绘图的函数写成一个完整的字符串。 >>f=‘ 2*exp(-x) .* sin(x) ‘; 式中，运用数组乘法定义了函数 f x e x x ( ) = sin( ) − 2 >>fplot(f , [0 8]); >>title(f) , xlabel(‘x‘)
2'exp(-×).*sin(x) 0.71 0.2 图132f(x)=2e-sin(x)的曲线 在区间0≤x≤8绘出上述函数,产生如图13.2所示的图形 除了这些基本特性,函数 fplot还有很多强大的功能,有关详细的信息,参阅《 MATLAB 参考指南》或在线帮助。 132极小化 作图除了提供视觉信息外,还常常需要确定一个函数的其它更多的特殊属性。在许多 应用中,特别感兴趣的是确定函数的极值,即最大值(峰值)和最小值(谷值)。数学上, 可通过确定函数导数(斜率)为零的点,解析上求出这些极值点。检验 humps的图形在峰 值和谷值点上的斜率就很容易理解这个事实。显然,如果定义的函数简单,则这种方法常常 奏效。然而,即使很多容易求导的函数,也常常很难找到导数为零的点。在这种情况下,以 及很难或不可能解析上求得导数的情况下,必须数值上寻找函数的极值点。 MATLAB提供 了两个完成此功能的函数fmin和 fmins。这两个函数分别寻找一维或n维函数的最小值。 这里仅讨论fmin。有关 fmins的详细信息,参阅《 MATLAB参考指南》。因为f(x)的最大 值等于-fx)的最小值,所以,上述fmin和 fmins可用来求最大值和最小值。如果还不清楚, 把上述图形倒过来看,在这个状态下,峰值变成了谷值,而谷值则变成了峰值。 为了解释求解一维函数的最小值和最大值,再考虑上述例子。从图13.2可知,在 附近有一个最大值,并且在xmm=4附近有一个最小值。而这些点的解析值为 Xma=兀/4≈0.785和Xmin=5m/4≈393。为了方便,用文本编辑器编写一个脚本M文 件,并用fmin寻出数值上极值点,给出函数主体如下 2*exp(-X)*sin(x) %o define function for min xmin=fmin(fn, 2, 5) lo search over range 2<x<5
0 2 4 6 8 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 x 2*exp(-x).*sin(x) 图 13.2 f x e x x ( ) = sin( ) − 2 的曲线 在区间 0 x 8 绘出上述函数，产生如图 13.2 所示的图形。 除了这些基本特性，函数 fplot 还有很多强大的功能，有关详细的信息，参阅《MATLAB 参考指南》或在线帮助。 13.2 极小化 作图除了提供视觉信息外，还常常需要确定一个函数的其它更多的特殊属性。在许多 应用中，特别感兴趣的是确定函数的极值，即最大值（峰值）和最小值（谷值）。数学上， 可通过确定函数导数（斜率）为零的点，解析上求出这些极值点。检验 humps 的图形在峰 值和谷值点上的斜率就很容易理解这个事实。显然，如果定义的函数简单，则这种方法常常 奏效。然而，即使很多容易求导的函数，也常常很难找到导数为零的点。在这种情况下，以 及很难或不可能解析上求得导数的情况下，必须数值上寻找函数的极值点。MATLAB 提供 了两个完成此功能的函数 fmin 和 fmins。这两个函数分别寻找一维或 n 维函数的最小值。 这里仅讨论 fmin。有关 fmins 的详细信息，参阅《MATLAB 参考指南》。因为 f(x)的最大 值等于-f(x)的最小值，所以，上述 fmin 和 fmins 可用来求最大值和最小值。如果还不清楚， 把上述图形倒过来看，在这个状态下，峰值变成了谷值，而谷值则变成了峰值。 为了解释求解一维函数的最小值和最大值，再考虑上述例子。从图 13.2 可知，在 xmax=0.7 附近有一个 最大值， 并且在 xmin=4 附近 有一个最 小值。而 这些点的 解析值为： xmax = / 4 0.785 和 xmin = 5 / 4 3.93 。为了方便，用文本编辑器编写一个脚本 M 文 件，并用 fmin 寻出数值上极值点，给出函数主体如下： % ex_fmin.m fn=‘ 2*exp(-x)*sin(x) ‘; % define function for min xmin=fmin(fn , 2 , 5) % search over range 2<x<5
emin=5*pi/4-xmin lo find error x-Xmin, need x since fn has x as its variable ymin=eval(fn) evaluate at xmin fx=-2*exp(-x)*sin(x) define for max note minus sign maximin( fx, 0, 3) lo search over range 0>ex-fmin Xmin- 3.9270 emin 14523e-006 0.0279 xmax- 0.7854 1.378le-005 ymax- 0.6448 这些结果与上述图形非常吻合。注意,fmin的工作方式很像 fplot。要计算的函数可用 函数M文件表达,或者只给出一个x为自变量的字符串。上述例子就是使用后一种方 法。这个例子也使用了函数eval,它获取一个字符串,并解释它,如同在 MATLAB提示符 下输入该字符串。由于要计算的函数以x为自变量的字符串形式给出,那么设置x等于xmin 和xmax,允许eva计算该函数,找到ymin和ymax。 最后,特别注意,求数值上的最小值包含一个搜索过程,fmin不断计算函数值,寻求 其最小值。如果计算的函数需要很大的计算量,或者该函数在搜索区间不止一个最小值,则 该搜索过程所花的时间比较长。在有些情况下,搜索过程根本找不到结果。当fmin找不到 最小值时,它会停止运行并提供解释 与函数fmin一样,函数 fmins搜索最小值。不过, fmins搜索向量的标量函数的最小 值。即 fmins寻找 f(x) 这里x是函数()的向量参数,函数f()返回标量值。函数 fmins利用单纯形法求最小值,它 不需要精确的梯度计算。任何一种优化工具箱中具有更多扩展的优化算法
emin=5*pi / 4-xmin % find error x=xmin; % need x since fn has x as its variable ymin=eval(fn) % evaluate at xmin fx=‘ -2*exp(-x)*sin(x) ‘; % define for max:note minus sign xmax=fmin(fx , 0 , 3) % search over range 0>ex-fmin xmin = 3.9270 emin = 1.4523e-006 ymin = -0.0279 xmax = 0.7854 emax = -1.3781e-005 ymax = 0.6448 这些结果与上述图形非常吻合。注意，fmin 的工作方式很像 fplot。要计算的函数可用 一个函数 M 文件表达，或者只给出一个 x 为自变量的字符串。上述例子就是使用后一种方 法。这个例子也使用了函数 eval，它获取一个字符串，并解释它，如同在 MATLAB 提示符 下输入该字符串。由于要计算的函数以 x 为自变量的字符串形式给出，那么设置 x 等于 xmin 和 xmax，允许 eval 计算该函数，找到 ymin 和 ymax。 最后，特别注意，求数值上的最小值包含一个搜索过程，fmin 不断计算函数值，寻求 其最小值。如果计算的函数需要很大的计算量，或者该函数在搜索区间不止一个最小值，则 该搜索过程所花的时间比较长。在有些情况下，搜索过程根本找不到结果。当 fmin 找不到 最小值时，它会停止运行并提供解释。 与函数 fmin 一样，函数 fmins 搜索最小值。不过，fmins 搜索向量的标量函数的最小 值。即 fmins 寻找 x min f(x) 这里 x 是函数 f(.)的向量参数，函数 f(.)返回标量值。函数 fmins 利用单纯形法求最小值，它 不需要精确的梯度计算。任何一种优化工具箱中具有更多扩展的优化算法
13.3求零点 正如人们对寻找函数的极点感兴趣一样,有时寻找函数过零或等于其它常数的点也非 常重要。一般试图用解析的方法寻找这类点非常困难,而且很多时候是不可能的。在上述函 数 humps的图中(如图13.3所示),该函数在x=12附近过零。 FPLOT OF HUMPS 0.5 1.5 图13.3 humps函数的图形 MATLAB再一次提供了该问题的数值解法。函数 fzero寻找一维函数的零点。为了说 明该函数的使用,让我们再运用 humps例子 >>xzero-fzero(humps, 1.2) look for a zero near 1.2 Zero- >>zero=humps(zero, 1.2) lo evaluate at zero 3.5527e-15 所以, humps的零点接近于1.3。如前所述,寻找零点的过程可能失败。如果 fzero没有找到 零点,它将停止运行并提供解释。 当调用函数 fzero时,必须给出该函数的名称。但由于某种原因,它不能接受以x为自 变量的字符串来描述的函数。这样,即使在 fplot和fmin中都具有的这个特性, fzero将不 fzero不仅能寻找零点,它还可以寻找函数等于任何常数值的点。仅仅要求一个简单的
13.3 求零点 正如人们对寻找函数的极点感兴趣一样，有时寻找函数过零或等于其它常数的点也非 常重要。一般试图用解析的方法寻找这类点非常困难，而且很多时候是不可能的。在上述函 数 humps 的图中(如图 13.3 所示)，该函数在 x=1.2 附近过零。 0 0.5 1 1.5 2 -20 0 20 40 60 80 100 FPLOT OF HUMPS 图 13.3 humps 函数的图形 MATLAB 再一次提供了该问题的数值解法。函数 fzero 寻找一维函数的零点。为了说 明该函数的使用，让我们再运用 humps 例子。 >>xzero=fzero(‘ humps ‘ , 1.2) % look for a zero near 1.2 xzero= 1.2995 >>yzero=humps(xzero , 1.2) % evaluate at xzero yzero= 3.5527e-15 所以，humps 的零点接近于 1.3。如前所述，寻找零点的过程可能失败。如果 fzero 没有找到 零点，它将停止运行并提供解释。 当调用函数 fzero 时，必须给出该函数的名称。但由于某种原因，它不能接受以 x 为自 变量的字符串来描述的函数。这样，即使在 fplot 和 fmin 中都具有的这个特性，fzero 将不 工作。 fzero 不仅能寻找零点，它还可以寻找函数等于任何常数值的点。仅仅要求一个简单的
再定义。例如,为了寻找f(x=c的点,定义函数g(x)=f(x)c,然后,在fero中使用gx), 就会找出g(x)为零的x值,它发生在f(x)=c时。 13.4积分 个函数的积分或面积也是它的另一个有用的属性。 MATLAT提供了在有限区间内 数值计算某函数下的面积的三种函数:trap2,quad和quad8。函数trap通过计算若干梯 形面积的和来近似某函数的积分,这些梯形如图134所示,是通过使用函数 humps的数据 点形成 彐40 15005 15 图134粗略的梯形逼近曲线下的面积示意图 从图中可明显地看出,单个梯形的面积在某一段欠估计了函数真正的面积,而在其它 段又过估计了函数的真正面积。如同线性插值,当梯形数目越多时,函数的近似面积越准确 例如,在图134中,如果我们大致增加一倍数目的梯形,我们得到如下页(如图135)所 示的更好的近似结果
再定义。例如，为了寻找 f(x)=c 的点，定义函数 g(x)=f(x)-c，然后，在 fzero 中使用 g(x)， 就会找出 g(x)为零的 x 值，它发生在 f(x)=c 时。 13.4 积分 一个函数的积分或面积也是它的另一个有用的属性。MATLAT 提供了在有限区间内， 数值计算某函数下的面积的三种函数：trap2 , quad 和 quad8。函数 trapz 通过计算若干梯 形面积的和来近似某函数的积分，这些梯形如图 13.4 所示，是通过使用函数 humps 的数据 点形成。 图 13.4 粗略的梯形逼近曲线下的面积示意图 从图中可明显地看出，单个梯形的面积在某一段欠估计了函数真正的面积，而在其它 段又过估计了函数的真正面积。如同线性插值，当梯形数目越多时，函数的近似面积越准确。 例如，在图 13.4 中，如果我们大致增加一倍数目的梯形，我们得到如下页（如图 13.5）所 示的更好的近似结果
60 0.5 15 图13.5较好的梯形逼近曲线下的面积示意图 对如上所示的两个曲线,用 trapz在区间-1>X=-1:0.17:2; rough approximated >>area=trapz(x, y) %o call trapz just like the plot command 25.9174 >>X=-1:0.07:2; %o better approximate >>area=trapz(x, y) 26.6243 自然地,上述两个结果不同。基于对图形的观察,粗略近似可能低估了实际面积。除 非特别精确,没有准则说明哪种近似效果更好。很明显,如果人们能够以某种方式改变单个 梯形的宽度,以适应函数的特性,即当函数变化快时,使得梯形的宽度变窄,这样就能够得 到更精确的结果 MATLAB的函数quad和quad8是基于数学上的正方形概念来计算函数的面积,这些 积分函数的操作方式一样。为获得更准确的结果,两个函数在所需的区间都要计算被积函数 此外,与简单的梯形比较,这两个函数进行更高阶的近似,而且quad8比quad更精确。这 两个函数的调用方法与fero相同, >>area=quad( humps,-1, 2)% find area between-I and 2
图 13.5 较好的梯形逼近曲线下的面积示意图 对如上所示的两个曲线，用 trapz 在区间-1>x=-1 : 0.17 : 2; % rough approximation >>y=humps(x); >>area=trapz(x , y) % call trapz just like the plot command area = 25.9174 >>x=-1 : 0.07 : 2; % better approximation >>y=humps(x); >>area=trapz(x , y) area = 26.6243 自然地，上述两个结果不同。基于对图形的观察，粗略近似可能低估了实际面积。除 非特别精确，没有准则说明哪种近似效果更好。很明显，如果人们能够以某种方式改变单个 梯形的宽度，以适应函数的特性，即当函数变化快时，使得梯形的宽度变窄，这样就能够得 到更精确的结果。 MATLAB 的函数 quad 和 quad8 是基于数学上的正方形概念来计算函数的面积，这些 积分函数的操作方式一样。为获得更准确的结果，两个函数在所需的区间都要计算被积函数。 此外，与简单的梯形比较，这两个函数进行更高阶的近似，而且 quad8 比 quad 更精确。这 两个函数的调用方法与 fzero 相同，即 >>area=quad(‘ humps ‘ , -1 , 2) % find area between -1 and 2
irea- >>areaquad8(humps",-1, 2) area 26.3450 注意,这两个函数返回完全相同的估计面积,而且这个估计值在两个 trapz面积的估计 值之间。有关 MATLAB的积分函数的其它信息,参阅《 MATLAB参考指南》或在线帮助。 135微分 与积分相反,数值微分非常困难。积分描述了一个函数的整体或宏观性质,而微分则 描述一个函数在一点处的斜率,这是函数的微观性质。因此积分对函数的形状在小范围内的 改变不敏感。而微分却很敏感。一个函数小的变化,容易产生相邻点的斜率的大的改变。 由于微分这个固有的困难,所以尽可能避免数值微分,特别是对实验获得的数据进行 微分。在这种情况下,最好用最小二乘曲线拟合这种数据,然后对所得到的多项式进行微分。 或用另一种方法,对该数据进行三次样条拟合,然后寻找如第11章所讨论的样条微分。例 如,再次考虑第11章曲线拟合的例子。 >y=-.4471.9783.286.167.087.347669.569.489.3011.2];%data >>n=2: order of fit >>p=polyfit(x, y, n) find polynomial coefficients 9.810820.1293-0.0317 >>Zpolyval(p, xi); evaluate polynomial >>xlabel(x), ylabel( y=f(x)), title( Second Order Curve Fitting) 在这种情况下,运用多项式微分函数 polder求得微分。 19621720.1293
area = 26.3450 >>area=quad8(‘ humps ‘ , -1 , 2) area = 26.3450 注意，这两个函数返回完全相同的估计面积，而且这个估计值在两个 trapz 面积的估计 值之间。有关 MATLAB 的积分函数的其它信息，参阅《MATLAB 参考指南》或在线帮助。 13.5 微分 与积分相反，数值微分非常困难。积分描述了一个函数的整体或宏观性质，而微分则 描述一个函数在一点处的斜率，这是函数的微观性质。因此积分对函数的形状在小范围内的 改变不敏感。而微分却很敏感。一个函数小的变化，容易产生相邻点的斜率的大的改变。 由于微分这个固有的困难，所以尽可能避免数值微分，特别是对实验获得的数据进行 微分。在这种情况下，最好用最小二乘曲线拟合这种数据，然后对所得到的多项式进行微分。 或用另一种方法，对该数据进行三次样条拟合，然后寻找如第 11 章所讨论的样条微分。例 如，再次考虑第 11 章曲线拟合的例子。 >>x=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1] >>y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; % data >>n=2; % order of fit >>p=polyfit(x , y , n) % find polynomial coefficients p = -9.8108 20.1293 -0.0317 >>xi=linspace(0 , 1 , 100); >>z=polyval(p , xi); % evaluate polynomial >>plot(x , y , ‘ o ' , x , y , xi , z , ' : ') >>xlabel(‘ x ‘) , ylabel(‘ y=f(x) ‘) , title(‘ Second Order Curve Fitting ‘) 在这种情况下，运用多项式微分函数 polyder 求得微分。 >>pd=polyder(p) pd = -19.6217 20.1293
Second Order Curve Fitting 10 -2 0.2 0.6 08 图136二次曲线拟合 y=-98108x2+201293x-00317的微分是dydx=19.6217x+2012093由于一个多 项式的微分是另一个低一阶的多项式,所以还可以计算并画出该函数的微分。 >>z=polyval(pd, xi); evaluate derivative >>plot(xi, z) >>xlabel( x), ylabel( dy/dx), title( Derivative of a curve Fit Polynimial) (微分曲线如图13.7所示)
图 13.6 二次曲线拟合 y = −9 8108x + x − 2 . 20.1293 0.0317 的微分是 dy/dx=-19.6217x+20.1293。由于一个多 项式的微分是另一个低一阶的多项式，所以还可以计算并画出该函数的微分。 >>z=polyval(pd , xi); % evaluate derivative >>plot(xi , z) >>xlabel(‘ x ‘) , ylabel(‘ dy/dx ‘) , title(‘ Derivative of a curve Fit Polynimial ‘) (微分曲线如图 13.7 所示）
Derivative of a curve Fit polynimial 谷 0 0 0.2 04 0.6 0.8 图137曲线拟合多项式微分 在这种情况下,拟合的多项式为二阶,使其微分为一阶多项式。这样,微分为一条直 线,它意味该微分与x成线性变化。 给定一些描述某函数的数据, MATLAB提供了一个计算其非常粗略的微分的函数。这 个函数命名为dif,它计算数组中元素间的差分。因为微分定义为: dy f(x+h)-f(x) dxh→0(x+h)-(x) 则y=f(x)的微分可近似为: dy f(x+h)-f(x) 这里h>0 dx (x+h)-(x) 它是y的有限差分除以x的有限差分。因为diff计算数组元素间的差分,所以在 MATLAB 中,可近似求得函数的微分。继续前一个例子 dy=diff(y). /diff(x); compute differences and use array divisie >>xd=x(1: length(x)-1); create new x axis since dy is shorter than y >>plot(xd, dy); >>title( Approximate Derivative Using DIFF >>ylabel( dy/dx), xlabel(x)
0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 x dy/dx Derivative of a curve Fit Polynimial 图 13.7 曲线拟合多项式微分 在这种情况下，拟合的多项式为二阶，使其微分为一阶多项式。这样，微分为一条直 线，它意味该微分与 x 成线性变化。 给定一些描述某函数的数据，MATLAB 提供了一个计算其非常粗略的微分的函数。这 个函数命名为 diff，它计算数组中元素间的差分。因为微分定义为： dy dx f x h f x h x h x = + − →0 + − lim ( ) ( ) ( ) ( ) 则 y=f(x)的微分可近似为： dy dx f x h f x x h x + − + − ( ) ( ) ( ) ( ) 这里 h>0 它是 y 的有限差分除以 x 的有限差分。因为 diff 计算数组元素间的差分，所以在 MATLAB 中，可近似求得函数的微分。继续前一个例子： >>dy=diff(y) ./ diff(x); % compute differences and use array division >>xd=x(1 : length(x)-1); % create new x axis since dy is shorter than y >>plot(xd , dy); >>title(‘ Approximate Derivative Using DIFF ‘) >>ylabel(‘ dy/dx ‘) , xlabel(‘ x ‘)
