1、- 1 -第 7章 计算方法的 MATLAB实现计算方法主要研究数学问题的数值解，涉及的内容广。利用 MATLAB 提供的部分函数解决，可以实现某些情况下的数值求解。本章作简要介绍。7.1 一元非线性方程求解求解一元非线性方程的方法主要有二分法、割线法、牛顿法等。本节主要介绍两个可以求解一元非线性方程的函数，即 fzero 函数和 roots 函数。7.1.1 fzero 函数用 fzero 函数求一元非线性方程的零点。该函数的简单调用格式为：x=fzero(fun,x0)：如果 x0 为标量，则试图寻找函数 fun 在 x0 附近的零点。fun 参数是一个 M 文件函数或匿名函数的函数句柄。
2、函数 fzero 返回的值x 靠近 fun 函数改变符号的位置。如果搜索失败，则返回 NaN。搜索区间扩展到发现 inf、NaN 或复数值时终止搜索。如果 x0 是一个长度为 2 的矢量，则fzero 函数假设 x0 是一个区间，其中 fun(x0(1)的符号与 fun(x0(2)的不同。如果符号相同，则会出错。给出符号不同的区间可以保证 fzero 函数返回一个fun 函数改变符号的位置附近的值。x,fval=fzero()：返回解 x 和解 x 处目标函数 fun 的值。x,fval,exitflag=fzero()：还返回一个 exitflag 值，描述 fzero 函数的退出条件。返回
3、值及其描述如下所示：0，表示函数找到了零值点 x；0，表示没有发现零值点。例 71 求方程 的解。0)3(5)1(20程序：f=(x)(x+20)*(x+1)*(x-5)*(x-30);%匿名函数x1=fzero(f,-50)%解一元非线性方程x2=fzero(f,1)x3=fzero(f,6)x4=fzero(f,50)运行结果：x1 =-20x2 =- 2 -1x3 =5x4 =307.1.2 roots 函数用 roots 函数计算多项式的根，该函数的调用格式为：r=roots(c)：返回一个列矢量，其元素为多项式 c 的解。例 72 求多项式方程 的解。0523x程序：P=1 -5 1
4、 -5;%表示多项式 r=roots(P)%求多项式方程的解P1=poly(r)%由多项式的根返回多项式的系数运行结果：r =5.0000 -0.0000 + 1.0000i-0.0000 - 1.0000iP1 =1.0000 -5.0000 1.0000 -5.0000注意：poly 和 roots 互为逆函数，函数 poly 返回多项式的系数。7.2 线性方程组的数值解法用 MATLAB 求解线性方程组的解，在 6.3 中已介绍了基于矩阵变换的直接解法，除此方法之外，还有很多方法，比如，Jocabi 迭代法、Gauss-Seidel迭代法和 SOR(超松驰)迭代法等，这里不作详细介绍。可
5、参阅有关文献。7.3 非线性方程组的数值解法实际工程中得到的数学模型往往具有非线性的特点，得到其解析解比较困难。一般采用迭代法求解非线性方程组。比较常见的迭代方法有不动点迭代法、Newton 迭代法和拟 Newton 迭代法等几种。本节仅介绍不动点迭代法。设含有 个未知数和 个方程的非线性方程组记为： ，然后把n 0)(xF上述方程组改写为便于迭代的等价形式： ，由此得出不动点迭代法)(x- 3 -的迭代公式： 。)(1kkx如果得到的序列 满足 ，则 就是 的不动点。这样就*limxk可以求出线性方程组的解（近似解） 。据此，编写不动点迭代法的 M 文件（文件名为 example7_3b.m
6、）：function s=example7_3b(x,eps)%不动点迭代法求非线性方程组的解%x 为迭代初值，eps 为允许误差值。if nargin=1eps=1e-6;elseif nargin=epsx=x1;x1=example7_3a(x);%循环迭代ends=x1;return例 73 用不动点迭代法求方程组 的解。0869122xx解：变形为 /)8(91221x函数代码（M 文件名为 example7_3a.m）：function y=example7_3a(x)y(1)=(x(1)*x(1)+x(2)*x(2)+9)/(-11);y(2)=(x(1)*x(2)*x(2)+x
7、(1)+8)/6;调用函数函数格式：s=example7_3b(0,0)运行结果：s=-1.0000 1.00007.4 插值- 4 -插值计算在数据拟合和数据平滑等方面应用普遍。MATLAB 提供了用最近邻插值、线性插值、三次样条插值、三次插值和 FFT 插值法进行一维、二维、三维和高维插值的函数。7.4.1 一维插值MATLAB 中有两种一维插值，即多项式插值和基于 FFT 的插值。 多项式插值函数 interp1 进行一维插值。一维插值是进行数据分析和曲线拟合的重要手段。interp1 函数使用多项式技术，用多项式函数拟合所提供的数据，并计算目标插值点上的插值函数值，其最常用的语法形式是
8、：yi=interp1(x,y,xi,method)x 和 y 为给定数据的矢量，长度相同。xi 为包含要插值的点的矢量，method 是一个可选的字符串，指定一种插值方法，包括：最近邻插值（method=nearest）：该方法将插值点的值设置为已知数据点中距离最近的点的值；线性插值（method=linear ）：该方法用线性函数拟合每对数据点，并返回 xi 处的相关函数值；三次样条插值（method=spline）：该方法用三次样条函数拟合每对数据点，用 spline 函数在插值处进行三次样条插值；三次插值（method=pchip 或 cubic）：该方法用 pchip 函数对矢量 x
9、 和y 进行分段三次 Hermite 插值。这几种方法在速度、内存和平滑性等方面有所差别，使用时可以根据需要进行选择，包括：最近邻插值是最快的方法，但是，利用它得到的结果平滑性最差；线性插值比最近邻插值要占用更多的内存，运行时间略长。与最近邻法不同，它生成的结果是连续的，但是在顶点处有坡度变化；三次样条插值的运行时间相对来说最长，内存消耗比三次插值略少。它生成的结果平滑性最好。但是，如果输入数据不很均匀，可能会得到意想不到结果；三次插值需要更多内存，而且运行时间比最近邻插值和线性插值的长。但是，使用此法时，插值数据及其导数都是连续的。例 74程序：x=1:10;y=1800 778 518 5
10、000 980 588 2799 528 6700 598;x1=1:0.1:10;y1=interp1(x,y,x1,nearest);%一维最近邻插值- 5 -plot(x1,y1)y2=interp1(x,y,x1,linear);%一维线性插值hold onplot(x1,y2,color,r)hold off运行结果见图 71。图 71 一维最近邻插值和线性插值例 75程序：x=1:10;y=1800 778 518 5000 980 588 2799 528 6700 598;plot(x,y);x1=1:0.1:10;hold ony1=interp1(x,y,x1,spline
11、);%三次样条插值plot(x1,y1,color,r)hold off运行结果见图 72。- 6 -图 72 三次样条插值 基于 FFT 的插值函数 interpft 用基于 FFT 的方法进行一维插值。本方法计算包含周期函数值的矢量的傅里叶变换。然后，它用更多的点计算逆傅里叶变换。该函数的调用形式为：y=interpft(x,n)其中，x 是一个包含周期函数值的矢量，这些值在等间隔的点上采集。n 是样本大小。7.4.2 二维插值二维插值在图像处理和可视化方面有着很重要的应用。MATLAB 用函数interp2 进行二维插值。该函数的一般形式为：ZI=interp2(X,Y,Z,XI,YI,
12、method)其中，Z 是一个矩阵数组，包含二维函数的值，X 和 Y 为大小相同的数组，包含相对于 Z 的给定值。XI 和 YI 为包含插值点数据的矩阵， method 表示插值方法，为可选参数。MATLAB 提供了三种不同的插值方法进行二维插值：最近邻插值（method=nearest）：该方法用分区域常数曲面拟合数据，插值点的值是最近点的值；双线性插值（method=linear ）：该方法用双线性曲面拟合数据点，插- 7 -值点的值是四个最近的值的组合。本方法是分区域双线性的，比双三次插值法快，并且内存消耗更少；双三次插值（method=cubic）：该方法用双三次曲面拟合数据点，插值点
13、的值是 16 个最近点的值的组合。本方法是分区域三次的，结果的平滑性比前面两种的都好。注意：所有这些方法都要求 X 和 Y 数据是单调的，即从点到点，要么总是递增的，要么总是递减的。应该用 meshgrid 函数准备这些矩 阵。例 76程序：%低分辨率的 peaks 函数图形x,y=meshgrid(-4:1:4);z=peaks(x,y);surf(x,y,z)运行结果见图 73。图 73 低分辨率的 peaks 函数图形例 77程序：x,y=meshgrid(-4:1:4);z=peaks(x,y);xI,yI=meshgrid(-4:0.25:4);zI=interp2(x,y,z,xI
14、,yI,nearest);%二维最近邻插值surf(xI,yI,zI)运行结果见图 74。- 8 -图 74 二维最近邻插值例 78程序：x,y=meshgrid(-4:1:4);z=peaks(x,y);xI,yI=meshgrid(-4:0.25:4);zI=interp2(x,y,z,xI,yI,linear);%二维双线性插值,linear 可省略.surf(xI,yI,zI)运行结果见图 75。图 75 二维双线性插值- 9 -例 79程序：x,y=meshgrid(-4:1:4);z=peaks(x,y);xI,yI=meshgrid(-4:0.25:4);zI=interp2(x
15、,y,z,xI,yI,cubic);%二维双三次插值surf(xI,yI,zI)运行结果见图 76。图 76 二维双三次插值7.4.3 多维插值MATLAB 提供了几种多维数据的插值函数，如表 71 中所示。表 71 多维数据的插值函数函 数 描 述interp3 三维数据插值interpn 多维数据插值ndgrid 多维数据网格化函数 interp3 进行三维插值。计算三维样本集 V 中数据点之间的值。该函数的一般形式为：VI=interp3(X,Y,Z,V,XI,YI,ZI,method)其中，X，Y 和 Z 指定数据点；V 包含与 X，Y 和 Z 对应的值； XI，YI 和 ZI- 10
16、 -为 interp3 函数对 V 中数据进行插值的点。对于超出范围的值，interp3 函数返回 NaN。method 表示插值方法，为可选参数。对于三维数据，有三种不同的插值方法。最近邻插值（method=nearest）：该方法选择最近点的值；线性插值（method=linear ）：该方法基于最近的 8 个点进行分区域线性插值；三次插值（method=cubic ）：该方法基于最近的 64 个点进行分区域三次插值。用 interpn 函数进行更高维数据的插值，该函数的常用形式为：VI=interpn(X1,X2,X3,V,Y1,Y2,Y3,method)其中，Y1,Y2,Y3, 为 i
17、nterpn 函数对 V 中数据进行插值的点。对于超出范围的值，interpn 函数返回 NaN。method 表示插值方法，为可选参数。高维数据插值，同样有最近邻插值、线性插值和三次插值三种方法。用 ndgrid 函数为高维函数评价和插值生成数据数组。该函数将一系列输入矢量指定的图域转换为一系列输出数组。第 i 维是输入矢量 xi 的拷贝。ndgrid 函数的语法格式为：X1,X2,X3, = ndgrid(x1,x2,x3,)7.5 曲线拟合所谓曲线拟合，就是利用两个或多个变量的离散点，用平滑曲线来拟合它们之间的关系。根据拟合方法的不同，有参数拟合和非参数拟合之分。参数拟合，曲线不要求通过
18、所有点，采用最小二乘法；非参数拟合，要求曲线通过所有点，采用插值法。由于曲线拟合是数据分析最常见的任务之一，MATLAB提供了多种函数和工具来进行曲线拟合，另外还有曲线拟合工具箱。7.5.1 最小二乘法最小二乘法通过最小化残差的平方和来获得待定系数的估计。第 个数据i点的残差定义为测量响应值 和拟合响应值 之间的差值，即 ，iyiyiiyri残差的平方和 。niiinirS1212)(常见的最小二乘法包括线性最小二乘、加权线性最小二乘、稳健最小二乘和非线性最小二乘等。求解非线性最小二乘问题的 Gauss-Newton 法和Levenberg-Marquart 法是老牌算法。- 11 -7.5.
19、2 多项式曲线拟合用 polyfit 函数计算拟合数据集的多项式在最小二乘意义上的系数，调用形式为：p=polyfit(x,y,n)x 和 y 是包含要拟合的 和 数据的矢量，n 是多项式的阶次。xy例 710程序：x=1:10;y=1 3 31 133 381 871 1723 3081 5113 8011;P=polyfit(x,y,4)%多项式曲线拟合，返回多项式的系数。运行结果：P =1.0000 -2.0000 -0.0000 1.0000 1.00007.5.3 相关工具MATLAB 支持用基本拟合界面进行拟合。该拟合界面具有拟合快速，操作简便的优势。它具有如下功能：使用样条插值、
20、分段三次艾尔米特插值（PCHIP）或者是 1 到 10 阶的多项式插值进行数据拟合；利用一组数据可以同时作多条拟合曲线；可以绘制残差图；可以检查拟合结果；可以对拟合进行内插或外推；用拟合结果和标准残差在图中进行注释；可以将拟合和计算结果保存到 MATLAB 工作空间。下面结合一个具体的例子加以说明。例 711 某商店某种产品的销售量如表 72 所示。表 72 某产品销售量资料年份 2001 2002 2003 2004 2005 2006 2007 2008 2009销售量（万件） 10.0 18.0 25.0 30.5 35.0 38.0 40.0 39.5 38.0请用多项式曲线拟合上述数
21、据。按照下面的步骤进行操作。用上述数据绘散点图和线形图的组合图（见图 77） ；在图形窗口的“Tools”菜单中单击“Basic Fitting”菜单选项；- 12 -两次单击“”按钮。打开的曲线拟合界面“Basic Fitting”对话框如图 78 所示。基本拟合界面中各选项的功能包括：Select data 下拉式列表框：该列表框中显示了图形窗口中图形用到的所有数据集的名称，在其中选择要拟合的数据。一次只能选择一组数据，但在一组数据里可以同时拟合多条曲线。Center and scale x data 单选钮：选择此项以后，数据中心化为具有 0 均值，比例化为具有单位标准差。对数据进行中心
22、化和比例化，可以提高后面数值计算的精度。图 77 散点图和线形图的组合图Plot fits：使用本面板，可以用图形查看当前数据集的一种或多种拟合结果。Check to display fits on figure：选择当前数据集的拟合类型。有两种拟合类型可供选择，即插值和多项式。进行三次样条插值使用 Spline 函数，保形（shape-preserving）插值使用 pchip 函数（三次插值） 。多项式拟合使用polyfit 函数。可以选择多种拟合类型。Show equations 单选钮：选此项，在拟合图形上显示方程。- 13 -Plot residuals 单选钮：选此项，显示拟合曲线
23、的残差，可用条形图、散点图或线形图显示。Show norm of residuals 单选钮：选此项，显示残差的范数。残差的范数是表示拟合优度的一个统计量，值越小，表示拟合程度越高。用 norm 函数进行计算，即 norm(V,2)，其中 V 为残差矢量。图 78 基本拟合界面Numerical results 方框：使用该面板，可以在不绘拟合图的情况下探察对当前数据集进行单次拟合的数值结果。Fit：选择拟合当前数据集的方程。拟合结果显示在菜单下面的列表框中。注意，在该菜单中选择方程并不影响“Plot fits”面板的状态。所以。如果试图在数据散点图中显示拟合曲线，可能需要在“Plot fit
24、s”面板中选择有关的核选框。Coefficients and norm of residuals：显示“Fit”中选择的方程的计算结果。Save to workspace：打开一个对话框，使用它将拟合结果保存到工作空间变量中。Find Y=f(X)：对当前的拟合进行内插或外推。Enter value(s)：输入一个 MATLAB 表达式，进行拟合计算。单击- 14 -“Evaluate”按钮以后计算表达式，结果显示在有关的表格中。当前拟合显示在“Fit”菜单中。Save to workspace：打开一个对话框，使用它将拟合结果保存到工作空间变量中。Plot evaluated results
25、：选择此项，计算结果显示到图中的数据点上。在基本拟合界面中作以下设置（见图 79）：用二次多项式拟合数据；在拟合图上显示方程；将拟合残差作为条形图显示，并将残差的图形作为子图显示；显示残差的范数。当前数据集用二次多项式拟合数据显示方程绘残差图显示残差的范数图 79 进行选项设置利用“Plot fits”面板可以可视地探察当前数据集的多个拟合图形。为了进行比较，通过选择合适的核选框来拟合数据的其他方程。如果某个方程生成的结果在数值上不精确，MATLAB 会显示一个警告信息。此时，应该选择“Center and scale x data” 核选框来改进数值精度。最后将拟合结果和残差显示在图 710
26、 中。图例显示了数据集和方程的名称。如果图例覆盖了图形的一部分，可以通- 15 -过单击和拖拉操作将它移到其他地方。图 710 拟合结果和残差图可以指定一个包含 x 值的矢量，计算这些位置上的拟合值。将该矢量输入到“Evaluate ”按钮左边的文本框中，然后单击“Evaluate”按钮。例如，输入矢量 2010:2015 后，再单击 “Evaluate”按钮，即可得到 2010 年到 2015 年的销售量的估计。如图 711 所示。- 16 -图 711 显示 x 值和对应的拟合值选择“Plot evaluated results”核选框，显示当前图形中对应当前数据的计算值，如图 712 所示图 712 拟合结果和对应的残差图习题七1、完成实验指导书中的实验六。
