第9章 数学建模函数
知识要点
本书前面几章已经详细介绍了MATLAB的各种基础知识及编程，本章重点介绍常用于MATLAB建模的函数，包括曲线拟合函数、参数估计函数等内容。
学习要求
9.1 曲线拟合函数
在科学和工程领域，曲线拟合的主要功能是寻求平滑的曲线来最好地表现带有噪声的测量数据，从这些测量数据中寻求两个函数变量之间的关系或变化趋势，最后得到曲线拟合的函数表达式y = f(x)。
使用多项式进行数据拟合会出现数据振荡，而Spline插值的方法可以得到很好的平滑效果，但是关于该插值方法有太多的参数，不适合用作曲线拟合的方法。
同时，由于在进行曲线拟合的时候，已经认为所有测量数据中已经包含噪声，因此，最后的拟合曲线并不要求通过每一个已知数据点，衡量拟合数据的标准则是整体数据拟合的误差最小。
一般情况下，MATLAB的曲线拟合方法用的是“最小方差”函数，其中方差的数值是拟合曲线和已知数据之间的垂直距离。
9.1.1 多项式拟合
在MATLAB中，函数polyfit()采用最小二乘法对给定的数据进行多项式拟合，得到该多项式的系数。该函数的调用格式如下。
● polyfit(x,y,n)：找到次数为n的多项式系数，对于数据集合{(xi ,yi)}，满足差的平方和最小。
● [p,E]=polyfit(x,y,n)：返回同上的多项式p和矩阵E。多项式系数在向量p中，矩阵E用在polyval函数中来计算误差。
例9-1：某数据的横坐标为x=[0.2 0.3 0.5 0.6 0.8 0.9 1.2 1.3 1.5 1.8]，纵坐标为y=[1 2 3 5 6 7 6 5 4 1]，对该数据进行多项式拟合。
在命令行窗口中依次输入：
clear all x = [0.3 0.4 0.7 0.9 1.2 1.9 2.8 3.2 3.7 4.5]; y = [1 2 3 4 5 2 6 9 2 7]; p5 = polyfit(x, y, 5); % 5阶多项式拟合 y5 = polyval(p5, x); p5 = vpa(poly2sym(p5), 5) % 显示5阶多项式 p9=polyfit(x, y, 9); % 9阶多项式拟合 y9=polyval(p9, x); figure; % 画图显示 plot(x, y, 'bo'); hold on; plot (x, y5, 'r:'); plot (x, y9, 'g--'); legend('原始数据', '5阶多项式拟合' ,'9阶多项式拟合'); xlabel ('x'); ylabel ('y');
运行程序后，得到的5阶多项式如下：
p5 = 0.8877*x^5 - 10.3*x^4 + 42.942*x^3 - 77.932*x^2 + 59.833*x - 11.673 警告: 多项式未正确设置条件。请添加具有不同 X 值的点，减少多项式的次数，或者尝试按照 HELP POLYFIT 中所述进行中心化和缩放。 > In polyfit (line 72) In ex9_01 (line 7)
运行程序后，得到的输出结果如图9-1所示。由图9-1可以看出，在使用5阶多项式拟合时，得到的结果比较差。
图9-1 多项式曲线拟合结果
当采用9阶多项式拟合时，得到的结果与原始数据比较相符。当使用函数polyfit()进行拟合时，多项式的阶次最大不超过length(x)-1。
9.1.2 加权最小方差（WLS）拟合原理及示例
所谓加权最小方差，就是根据基础数据本身各自的准确度的不同，在拟合的时候给每个数据以不同的加权数值。这种方法比前面介绍的单纯最小方差方法要更加符合拟合的初衷。
对应N阶多项式的拟合公式，所需要求解的拟合系数需要求解线性方程组，其中线性方程组的系数矩阵和需要求解的拟合系数矩阵分别为：
使用加权最小方差方法求解得到拟合系数：
其对应的加权最小方差为表达式
例9-2：根据WLS数据拟合方法，自行编写使用WLS方法拟合数据的M函数，然后使用WLS方法进行数据拟合。
创建M文件并命名为polyfits.m，利用M文件编辑器在M文件中输入：
function [th, err, yi] = polyfits(x, y, N, xi, r) % x, y: 数据点拟合 % N：多项式拟合的系统 % r: 加权系数的逆矩阵 M = length(x); x = x(:); y = y(:); % 判断调用函数的格式 if nargin == 4 % 当调用函数的格式为（x, y, N, r）时 if length(xi) == M r = xi; xi = x; % 当调用函数的格式为（x, y, N, xi）时 else r = 1; end % 当调用函数的格式为（x, y, N）时 elseif nargin == 3 xi = x; r = 1; end % 求解系数矩阵 A(:, N + 1) = ones(M, 1); for n = N : -1 : 1 A(:,n) = A(:,n + 1) .* x; end if length(r) == M for m = 1 : M A(m, :) = A(m, :) / r(m); y(m) = y(m) / r(m); end end % 计算拟合系数 th = (A \ y)'; ye = polyval(th, x); err = norm(y - ye) / norm(y); yi = polyval(th, xi);
使用上面的程序代码，对基础数据进行LS（最小方差）多项式拟合。在命令行窗口中依次输入：
clear all x = [-3 : 1 : 3]'; y = [1.1650 0.0754 -0.6965 0.0591 0.6268 0.3516 1.6961]'; [x, i] = sort(x); y = y(i); xi = min(x) + [0 : 100] / 100 * (max(x) - min(x)); for i = 1 : 4 N = 2 * i - 1; [th, err, yi] = polyfits(x, y, N, xi); subplot(2, 2, i) plot(x, y, 'o') hold on plot(xi ,yi, '-') grid on end
得到的拟合结果如图9-2所示。
图9-2 使用LS方法求解的拟合结果
从上面的例子中可以看出，LS方法其实是WLS方法的一种特例，相当于将每个基础数据的准确度都设为1。但是，自行编写的M文件和默认的命令结果不同，请仔细比较。
9.1.3 非线性曲线拟合
非线性曲线拟合是已知输入向量xdata、输出向量ydata，并知道输入与输出的函数关系为ydata=F(x,xdata)，但不清楚系数向量x。进行曲线拟合即求x使得下式成立：
在MATLAB中，可以使用函数curvefit解决此类问题，其调用格式如下。
● x=lsqcurvefit(fun,x0,xdata,ydata)：x0为初始解向量；xdata、ydata为满足关系ydata=F(x,xdata)的数据。
● x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub)：lb、ub为解向量的下界和上界，lb≤x≤ub，若没有指定界，则lb=[ ]，ub=[ ]。
● x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)：options为指定的优化参数。
● [x,resnorm]=lsqcurvefit(…)：resnorm是在x处残差的平方和。
● [x,resnorm,residual]=lsqcurvefit(…)：residual为在x处的残差。
● [x,resnorm,residual,exitflag]=lsqcurvefit(…)：exitflag为终止迭代的条件。
● [x,resnorm,residual,exitflag,output]=lsqcurvefit(…)：output为输出的优化信息。
例9-3：已知输入向量xdata和输出向量ydata，且长度都是n，使用最小二乘非线性拟合函数：
根据题意可知，目标函数为
其中：
初始解向量定位x0=[0.3, 0.4 ,0.1]。
首先创建M文件（拟合函数文件）并命名为ex9_03.m，利用M文件编辑器在M文件中输入：
function F = ex9_03(x, xdata) F = x(1) * xdata.^2 + x(2) * sin(xdata) + x(3) * xdata.^3;
再编写函数拟合代码：
clear all xdata = [3.6 7.7 9.3 4.1 8.6 2.8 1.3 7.9 10.0 5.4]; ydata = [16.5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3]; x0 = [10, 10, 10]; [x, resnorm] = lsqcurvefit(@ex9_03, x0, xdata, ydata)
结果为：
x = 0.2269 0.3385 0.3022 resnorm = 6.2950
即函数在x=0.2269、x=0.3385、x=0.3022处残差的平方和均为6.2950。
