数学建模预备知识——插值与拟合

插值与拟合

在实际问题中，一个函数 $y=f(x)$ 往往是通过实验观测得到的，仅已知函数 $f(x)$ 在某区间 $[a,b]$ 上一系列点的值，当需要这些节点 $x_0,x_1,\cdots ,x_n$ 之间的点上的函数值的时候，常用的简单的满足一定条件的函数 $\phi (x)$ 替代 $f(x)$，插值法就是一种常用的方法，其插值函数 $\phi (x)$ 满足
$$\phi (x_i)=y_i,i=1,\cdots ,n$$
拟合也是已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义上他的这些点的总偏差最小

插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上是完全不同的。

插值方法

在平面上给定一组离散的点列，要求用一条曲线把这些点按照次序连接起来，称为插值

分段线性插值

简单来说，将没两个相邻的节点用直线连接起来，如此形成一条折线就是分段线性插值函数，记作 $I_n(x)$，满足 $I_n(x_i)=y_i$，且在每个小区间 $[x_i,x_{i+1}]$ 上是线性函数

$I_n(x)$ 可以表示成 $I_n(x)=\sum _{i=0}^n{y}_i{l}_i(x)$，其中
$$l_i=\{\begin{array}{l}\frac{x-x_{i-1}}{x_i-x_{i-1}}x\in [x_{i-1},x_i],i\ne 0\\ \frac{x-x_{i+1}}{x_i-x_{i+1}}x\in [x_i,x_{i+1}],i\ne n\\ 0\end{array}$$
要求 $I_n(x)$ 具有良好的收敛性，即对于 $x\in [a,b]$ 有
$$\underset{x\to \infty }{\lim }{I}_n(x)=f(x)$$
用 $I_n(x)$ 计算 $x$ 点的插值时，只用到 $x$ 左右两个节点，计算量与节点数 $n$ 无关，但是 $n$ 越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了。

拉格朗日插值多项式

拉格朗日插值多项式的基函数为
KaTeX parse error: \cr valid only within a tabular/array environment
$l_i(x)$ 是一个 $n$ 次多项式，且满足
$$l_i(x_j)=\{\begin{array}{ll}0,& j\ne i\\ 1,& j=i\end{array}$$
拉格朗日插值函数为
$$\begin{gathered} L_n(x)=\sum _{i=0}^n{y}_i{l}_i(x)=\sum _{i=0}^n{y}_i(\prod _{j=0 \\ j\ne i}^n\frac{x-x_j}{x_i-x_j}) \end{gathered}$$

样条插值

在许多工程技术中，提出的计算问题对插值函数的光滑性有较高的要求，要求曲线不仅要连续，还要有连续的曲率，这就产生了样条插值

• 样条插值的概念

  样条本来是工程设计中的使用的一种绘图工具，是富有弹性的细木条或细金属条。绘图员利用它把一些已知的点连接成一条光滑的曲线，并使连接点处具有连续的曲率。三次样条插值就是由此抽象出来的。

  数学上，将具有一定光滑性的分段多项式称为样条函数。具体说，给定区间 $[a,b]$ 的一个分划
  $$\Delta :a=x_0<x_1<\cdots <x_{n-1}<x_n=b$$
  如果函数 $S(x)$ 满足

  1. 在每个小区间 $[x_i,x_{i+1}]$ 上都是 $S(x)$ 的 $m$ 次多项式
  2. $S(x)$ 在 $[a,b]$ 上具有 $m-1$ 阶连续导数

  则称 $S(x)$ 为关于分划 $\Delta$ 的 $m$ 次样条函数，其图形为 $m$ 次样条曲线

• 三次样条插值

  利用样条函数进行插值，即插值函数为样条函数，称为样条插值。

  下面介绍三次样条插值，即已知函数 $y=f(x)$ 在区间 $[a,b]$ 上的 $n+1$ 个节点
  $$a=x_0<x_i<\cdots <x_{n-1}<x_n=b$$
  上的值 $y_i=f(x_i)$ ，求插值函数 $S(x)$ ，使得

  1. $S(x_i)=y_i$
  2. 在每个小区间 $[x_i,x_{i+1}]$ 上，$S(x)$ 都是三次多项式，记为 $S_i(x)$
  3. $S(x)$ 在 $[a,b]$ 上二阶连续可微

  这样的函数 $S(x)$ 称为 $f(x)$ 的三次样条插值函数

  由条件2，不妨记
  $$\begin{array}{rl}& S(x)=S_i(x),x\in [x_i,x_{i+1}],i=0,1,\cdots ,n-1\\ & S_i(x)=a_i{x}^3+b_i{x}^2+c_{i}x+d\end{array}$$
  式中共有 $4n$ 个待定参数

  由条件3，不妨记
  $$\{\begin{array}{l}{S}_i(x_{i+1})=S_{i+1}(x_{i+1})\\ S_i^{\prime }(x_{i+1})=S_{i+1}^{\prime }(x_{i+1})\\ S_i^{\prime \prime }(x_{i+1})=S_{i+1}^{\prime \prime }(x_{i+1})\end{array}$$
  结合条件1和条件3，共能够得到 $4n-2$ 个方程，为了求出条件2的待定参数，还需要由两个边界条件

  常用的三次样条函数的边界条件由3中类型

  1. $S^{\prime }(a)=y_0^{\prime },S^{\prime }(b)=y_n^{\prime }$。

     由这种边界条件建立的样条插值函数称为 $f(x)$ 的完备三次样条插值函数。

     特别地，当 $y_0^{\prime }=y_n^{\prime }=0$ 时，样条曲线在端点处呈现水平状态

     如果 $f^{\prime }(x)$ 不知道，我们可以要求 $S^{\prime }(x)$ 与 $f^{\prime }(x)$ 在端点处近似相等，这时以 $x_0,x_1,x_2,x_3$ 为节点做一个三次 Newton 插值多项式 $N_a(x)$，以 $x_n,x_{n-1},x_{n-2},x_{n-3}$ 做一个三次 Newton 插值多项式 $N_b(x)$ ，要求
     $$S^{\prime }(a)=N^{\prime }(a),S^{\prime }(b)=N^{\prime }(b)$$
     由这种边界条件建立的三次样条称为 $f(x)$ 的Lagrange三次样条插值函数

  2. $S^{\prime \prime }(a)=y_0^{\prime \prime },S^{\prime \prime }(b)=y_n^{\prime \prime }$

     特别地，当 $y_0^{\prime \prime }=y_n^{\prime \prime }$ 时，称为自然边界条件

  3. $S^{\prime }(a+0)=S^{\prime }(b-0),S^{\prime \prime }(a+0)=S^{\prime \prime }(b-0)$

     此条件称为周期条件

Matlab插值工具箱

1. 一维插值函数

   Matlab中由现成的一维插值函数interp1，其基本语法为
   $$y=interp1(x_0,y_0,x{,}^{\prime }metho{d}^{\prime })$$
   其中 method指定插值的方法，其默认为线性插值。其全部参数如下

   参数	方法
   nearest	最近项插值
   linear	线性插值
   spline	立方样条
   cubic	立方插值

   所有的插值方法都要求 $x_0$ 是单调的

   当 $x_0$ 为等距时可以用快速插值法，使用快速插值法的格式为 *nearest``*linear``*spline``*cubic

2. 三次样条插值

   在Matlab中数据点称为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非纽结条件。这种条件强迫第1个和第2个三次多项式的三阶导数相等。对最后一个和倒数第二个两个三次多项式也做同样的处理。

   在Matlab中的三次样条插值有一下函数
   $$\begin{array}{rl}& y=interp1(x_0,y_0,x{,}^{\prime }splin{e}^{\prime })\\ & y=spline(x_0,y_0,x)\\ & pp=csape(x_0,y_0,x)orpp=csape(x_0,y_0,conds,valconds)\\ & y=fnval(pp,x)\end{array}$$
   其中：$x_0,y_0$ 为已知的数据点，$x,y$ 为插值点及其对应的函数值。

   对于三次样条插值，提倡使用函数csape。其返回值时pp形式，需要用fnval函数求得插值点的函数值。

   pp=scape(x0,y0)  
   

   使用的是默认的边界条件，即Lagrange边界条件

   pp=csape(x0,y0,conds,valconds)
   

   由conds指定插值的边界条件，其值有

   • complete 边界为一阶导数，一阶导数的值在valconds参数中给出
   • not-a-kont 非扭结条件
   • periodic 周期条件
   • second 边界为二阶导数，二阶导数的值在valconds参数中给出，如果忽略默认为 $[0.0]$
   • variational 设置边界的二阶导数值为 $[0,0]$

   对于一些特殊的边界条件，可以通过conds的一个 $1\times 2$ 矩阵来表示，conds 元素的取值为0，1，2。

   conds(i)=j的含义是给定端点 $i$ 的 $j$ 阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件，例：conds=[2,1]表示左边界为2阶导数，右边界为1阶导数，其值由valconds参数给出

   例：待加工零件的外形根据工艺要求由一组数据 $(x,y)$ 给出，加工时每一刀只能沿 $x$ 轴或者 $y$ 轴方向走非常小一步，这就需要从已知的数据得到加工所要求的步长的很小的 $(x,y)$ 坐标

   表中给出的 $x,y$ 数据位于机翼断面的下轮廓线上，假设需要得到 $x$ 坐标的每改变0.1时的 $y$ 坐标。试完成加工所需要的数据，画出曲线，并求出 $x=0$ 处的曲线斜率和 $13\le x\le 15$ 范围内 $y$ 的最小值。要求分别采用分段线性和三次样条两种插值方法

   $x$	0	3	5	7	9	11	12	13	14	15
   $y$	0	1.2	1.7	2.0	2.1	2.0	1.8	1.2	1.0	1.6

   解：

   x0=[0 3 5 7 9 11 12 13 14 15];
   y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
   x=0:0,1:15;
   y1=interp1(x0,y0,x,'linear');
   y2=interp1(x0,y0,x,'spline');
   pp1=csape(x0,y0);
   y3=fnval(pp1,x);
   pp2=csape(x0,y0,'second');
   y4=fnval(pp2,x);
   subplot(1,3,1)
   plot(x0,y0,'+',x,y1)
   title('Linear')
   subplot(1,3,2)
   plot(x0,y0,'+',x,y2)
   title('Spline1')
   subplot(1,3,3)
   plot(x0,y0,'+',x,y3)
   title('Spline2')
   ytmp=y3(130:150);
   ymin=min(ytmp);
   xmin=x(find(y3==ymin));
   [xmin,ymin]
   

   结果如图

   

3. 二维插值

   前面讲述的都是一维插值，即节点为一维变量，插值函数时一元函数。如果节点为二维，插值函数就是二元函数。如果在某区域测量若干节点的高程，为了画出较为精确的等高线图，要先插入更多的点，计算这些点的高程，这就是二维插值。

   1）插值节点为网格节点

   已知 $m\times n$ 个节点：$(x_i,y_j,z_{ij})$ ，且满足 $x_1<\cdots <x_m;y_i<\cdots <y_n$ 。求在 $(x,y)$ 处的插值 $z$。

   Matlab中有一些计算二维插值的命令
   $$z=interp2(x0,y0,z0,x,y{,}^{\prime }metho{d}^{\prime })$$
   其中：$x0,y0,z0$ 分别为 $m$ 维和 $n$ 维向量，表示节点；$z0$ 为 $n\times m$ 矩阵，表示节点值；$x,y$ 为一维数组，表示插值点，$x$ 与 $y$ 应是方向不同的向量，即一个是行向量一个是列向量。$z$ 是矩阵，行数为 $y$ 的维数，列数为 $x$ 的维数，表示得到的插值。

   如果是三次样条插值，可以使用命令
   p p = c s a p e ( { x 0 , y 0 } , z 0 , c o n d s , v a l c o n d s ) z = f n v a l ( p p , { x , y } ) \begin{aligned} &pp=csape(\{x0,y0\},z0,conds,valconds) \\ &z=fnval(pp,\{x,y\}) \end{aligned} ​pp=csape({x0,y0},z0,conds,valconds)z=fnval(pp,{x,y})​

   例：在一个丘陵地带测量高程，$x$ 和 $y$ 方向每隔100$m$ 测量一个点，得到高程如表所示，试插值求出一个曲面，确定合适的模型，并由此找出最高点和该点的高程

   

   解：

   x0=100:100:500;
   y0=100:100:400;
   z0=[636 697 624 478 450
   698 712 630 478 420
   680 674 598 412 400
   662 626 552 334 310];
   pp=csape({x0,y0},z0');
   x=100:1:500;
   y=100:1:400;
   z=fnval(pp,{x,y});
   zmax=max(max(z));
   [i,j]=find(z==max(max(z)));
   xx=x(i)
   yy=y(j)
   mesh(x,y,z')
   

   所得结果为：最高点是 $(166,177)$，高程为 $z=720.7540$。

   

   2）插值节点为散乱节点

   已知 $n$ 个节点 $(x_i,y_i,z_i)$ ，求点 $(x,y)$ 处的插值 $z$ 。

   对于上述问题，Matlab中提供了插值函数griddata
   $$ZI=griddata(x,y,z,XI,YI{,}^{\prime }metho{d}^{\prime })$$
   其中， $x,y,z$ 均为 $n$ 维向量，指出所给数据点的横坐标、纵坐标和竖坐标。向量 $XI,YI$ 是给定网格点的横坐标和纵坐标。返回值 $ZI$ 为网格 $(XI,YI)$ 处的函数值。 $XI$ 为行向量， $YI$ 为列向量

   参数	方法
   cubic	立方插值
   nearest	最近点插值

   P.S. Matlab插值时外插值有时是不确定的，可以采用立方插值和最近点插值混合插值法，把不确定的插值换成最近点插值。

   例：在某海域内测量得到一些点$(x,y)$ 处的水深 $z$ ，由下表给出，在适当的矩形区域内画出海底曲面的图形。

   

   解：

   x=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5];
   y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5];
   z=-[4 8 6 8 6 8 8 9 9 8 8 9 4 9];
   XI=min(x):max(x);
   YI=min(y):max(y);
   ZI1=griddata(x,y,z,XI,YI','cubic');
   ZI2=griddata(x,y,z,XI,YI','nearest');
   ZI=ZI1;
   ZI(isnan(ZI1))=ZI2(isnan(ZI1))
   mesh(XI,YI,ZI)
   

   

曲线拟合的线性最小二乘法

线性最小二乘法

线性拟合问题的提法是，已知一组数据，即平面上的 $n$ 个点 $(x_i,y_i)$ ，$x_i$ 互不相同，求一个函数 $y=f(x)$ ，使 $f(x)$ 在某种准则下与所有数据点最为接近，即曲线拟合效果最好

线性最小二乘法使解决曲线拟合的最常用的方法，其基本思路为：令
$$f(x)=a_1{r}_1(x)+a_2{r}_2(x)+\cdots +a_m{r}_m(x)$$
式中：$r_k(x)$ 为事先选定的一组线性无关的函数；$a_k$ 为待定系数

拟合准则是使 $y_i(i=1,2,\cdots ,n)$ 与 $f(x_i)$ 的距离 ${\delta }_i$ 的平方和最小，称为最小二乘准则。

1. 系数 $a_k$ 的确定

   记
   $$J(a_1,\cdots ,a_m)=\sum _{i=1}^n{\delta }_i^2=\sum _{i=1}^n[f(x_i)-y_i{]}^2$$
   为求 $a_1,\cdots ,a_m$ 使 $J$ 达到最小，只需要利用极值的必要条件 KaTeX parse error: Undefined control sequence: \part at position 7: \frac{\̲p̲a̲r̲t̲ ̲J}{\part a_j}=0… ，得到关于 $a_1,\cdots ,a_m$ 的线性方程组
   $$\sum _{i=1}^n{r}_j(x_i)[\sum _{k=1}^m{a}_k{r}_k(x_i)-y_i]=0,j=1,\cdots ,m$$
   即
   $$\sum _{k=1}^m{a}_k[\sum _{i=1}^n{r}_j(x_i)r_k(x_i)]=\sum _{i=1}^n{r}_j(x_i)y_{i}j=1,\cdots ,m$$
   记
   $$\begin{gathered} \mathbf{R}=\left[\begin{array}{ccc}{r}_1(x_1)& \cdots & r_m(x_1)\\ ⋮& & ⋮\\ r_1(x_n)& \cdots & r_m(x_n)\end{array}\right] \\ \mathbf{A}=[a_i,\cdots ,a_m{]}^T \\ \mathbf{Y}=[y_1,\cdots ,y_n{]}^T \end{gathered}$$
   由此上述方程组可以写成
   $${\mathbf{R}}^T\mathbf{R}\mathbf{A}={\mathbf{R}}^T\mathbf{Y}$$
   当 { r 1 ( x )   ⋯   , r m ( x ) } \{r_1(x)\,\cdots,r_m(x)\} {r1​(x)⋯,rm​(x)} 线性无关的时候，$\mathbf{R}$ 列满秩，${\mathbf{R}}^{\mathbf{T}}\mathbf{R}$ 可逆，于是方程组有唯一解
   $$\mathbf{A}=({\mathbf{R}}^{\mathbf{T}}\mathbf{R}{)}^{-1}{\mathbf{R}}^{\mathbf{T}}\mathbf{Y}$$

2. 函数 $r_k(x)$ 的选取

   面对一组数据 $(x_i,y_i)$ 用线性最小二乘法做曲线拟合的时候，最关键的一步就是恰当地选取 $r_1(x),\cdots ,r_m(x)$ 。如果通过机理分析，能够知道 $y$ 与 $x$ 之间的函数关系，则 $r_1(x),\cdots ,r_m(x)$ 容易确定。若无法知道 $y$ 与 $x$ 之间的函数关系，通常可以通过将数据 $(x_i,y_i)$ 作图，直观判断应该用什么曲线去拟合。通常有一下形式

   • 直线 $y=ax+b$
   • 多项式 $y=a_m{x}^m+\cdots +a_{1}x+a_0$
   • 双曲线 $y=\frac{a_1}{x}+a_2$
   • 指数曲线 $y=a_1{e}^{a_{2}x}$ （对于指数曲线，拟合前要取对数进行变量代换转化成线性函数）

最小二乘法的Matlab实现

1. 解方程组方法

   在上面的记号下，
   $$J(a_1,\cdots ,a_m)=\Vert \mathbf{R}\mathbf{A}-\mathbf{Y}{\Vert }_2^2$$
   Matlab中的线性最小二乘的标准型是
   $$\underset{A}{\min }\Vert \mathbf{R}\mathbf{A}-\mathbf{Y}{\Vert }_2^2$$
   系数矩阵 $\mathbf{A}=\mathbf{R}\backslash \mathbf{Y}$

   其中：$\mathbf{R}$ 为方程组的系数矩阵

   例：用最小二乘法求一个形如 $y=a+b{x}^2$ 的经验公式，使其与下表拟合

   x	19	25	31	38	44
   y	19.0	32.3	49.0	73.3	97.8

   解

   x0=[19 25 31 38 44]';
   y0=[19.0 32.3 49.0 73.3 97.8]';
   r=[ones(5,1),x0.^2];
   a=r\y0;
   x=19:0.1:44;
   y=a(1)+a(2)*x.^2;
   plot(x0,y0,'+',x,y,'r');
   

   所得经验公式为 $y=0.9726+0.05x^2$

   

2. 多项式拟合方法

   如果选取 { r 1 ( x ) , ⋯   , r m + 1 ( x ) } = { 1 , x , ⋯   , x m } \{r_1(x),\cdots,r_{m+1}(x)\}=\{1,x,\cdots,x^m\} {r1​(x),⋯,rm+1​(x)}={1,x,⋯,xm} ，即利用 $m$ 次多项式拟合给定数据，Matlab中现有的函数为polyfit ，具体命令如下
   $$p=polyfit(x0,y0,m)$$
   其中：$x0,y0$ 为要拟合的数据，$m$ 为多项式的次数，返回值 $p$ 为多项式的系数向量

   多项式在 $x$ 处的值可以用函数polyval 来计算，具体命令如下
   $$y=polyval(p,x)$$

   例：某企业近几年的生产利润如表所示

   年份	2010	2011	2012	2013	2014	2015	2016
   利润/万元	70	122	144	152	174	196	202

   试预测2017和2018年份的利润

   解：

   x0=2010:1:2016;
   y0=[70 122 144 152 174 196 202];
   plot(x0,y0,'*');
   p=polyfit(x0,y0,1);
   y2017=polyval(p,2017);
   y2018=polyval(p,2018);
   

   所得结果为 $y_{2017}=233.4286,y_{2018}=253.9286$

   

最小二乘优化

在无约束最优化问题中，有些重要的特殊情形，比如目标函数是由若干个函数的平方和构成的，这类函数一般可以写成
$$F(x)=\sum _{i=1}^m{f}_i^2(x),x\in R^n$$
式中：$x=[x_1,\cdots ,x_n{]}^T$

把极小化这类函数的问题
$$\min F(x)=\sum _{i=1}^m{f}_i^2(x)$$
称为最小二乘优化问题

最小二乘优化问题是一类比较特殊的优化问题，在处理这类问题的时候，Matlab也提供了一些强大的函数。在Matlab优化工具箱中，用于求解最小二乘优化问题的函数有lsqlin,lsqcurvefit,lsqnonlin,lsqnonneg

lsqlin函数

求解
$$\begin{gathered} \underset{x}{\min }\frac{1}{2}\Vert C\cdot x-d{\Vert }_2^2 \\ s.t.\{\begin{array}{l}A\cdot x\le b\\ A_{eq}\cdot x=b_{eq}\\ lb\le x\le ub\end{array} \end{gathered}$$
式中：$C,A,A_{eq}$ 均为矩阵，$d,b,b_{eq},lb,ub$ 为向量

在Matlab中的命令为
$$x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)$$

例：用最小二乘法求一个形如 $y=a+b{x}^2$ 的经验公式，使其与下表拟合

x	19	25	31	38	44
y	19.0	32.3	49.0	73.3	97.8

解：

x0=[19 25 31 38 44]';
y0=[19.0 32.3 49.0 73.3 97.8]';
r=[ones(5,1),x0.^2];
a=lsqlin(r,y);
x=19:0.1:44;
y=a(1)+a(2)*x.^2;
plot(x0,y0,'+',x,y,'r');

所得经验公式为 $y=0.9726+0.05x^2$

lsqcurvefit函数

给定输入输出数列 $xdata,ydata$ 求参量 $x$ ，使得
$$\underset{x}{\min }\Vert F(x,xdata)-ydata{\Vert }_2^2=\sum _i(F(x,xdat{a}_i)-ydat{a}_i{)}^2$$
Matlab中的函数为
$$x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)$$
其中：$fun$ 为定义函数 $F(x,xdata)$ 的.m文件

例：根据下列数据，拟合函数 $y=e^{-k_1{x}_1}sin(k_2{x}_2)+x_3^2$ 中的参数 $k_1,k_2$

$y/kg$	$x_1/c{m}^2$	$x_2/kg$	$x_3/kg$		$y/kg$	$x_1/c{m}^2$	$x_2/kg$	$x_3/kg$
15.02	23.73	5.49	1.21		15.94	23.52	5.18	1.98
12.62	22.34	4.32	1.35		14.33	21.86	4.86	1.59
14.86	28.84	5.04	1.92		15.11	28.95	5.18	1.37
13.98	27367	4.82	1.49		13.81	24.53	4.88	1.39
15.91	20.83	5.35	1.56		15.58	27.64	5.02	1.66
12.47	22.27	4.27	1.50		15.85	27.29	5.55	1.70
15.80	27.57	5.25	1.85		15.28	29.07	5.26	1.82
14.32	28.01	4.62	1.51		16.40	32.47	5.18	1.75
13.76	24.79	4.42	1.46		15.02	29.65	5.08	1.70
15.18	28.96	5.30	1.66		15.73	22.11	4.90	1.81
14.20	25.77	4.87	1.64		14.75	22.43	4.65	1.82
17.07	23.17	5.80	1.90		14.35	20.04	5.08	1.53
15.40	28.57	5.22	1.66

解：

先编写fun.m文件

function f=fun(x,xdata);
f=exp(-x(1)*xdata(:,1)).*sin(x(2)*xdata(:,2))+xdata(:,3).^2;

再编写主函数，将原始数据保存在文本文件data.txt中

a=textread('data.txt');
ydata=a(:,[2,7]);
ydata=nonzeros(ydata);
xdata=[a(:,[3:5]);a([1:end-1],[8:10])];
x=rand(2,1)
x=lsqcurvefit('fun',x,xdata,ydata,[],[]);

lsqnonlin函数

已知函数向量 $F(x)=[f_1(x),\cdots ,f_k(x){]}^T$ ，求 $x$ 使得
$$\underset{x}{\min }\Vert F(x){\Vert }_2^2$$
在Matlab中的函数为
$$x=lsqnonlin(fun,x0,lb,ub,options)$$
其中 $fun$ 为定义向量函数 $F(x)$ 的 .m 文件

例：用lsqnonlin函数拟合 $y=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$ 中的未知参数，其中的数据通过随机数生成。

x0=-10:0.1:10;
y0=normpdf(x0,0,1);
f=@(cs) 1/sqrt(2*pi)/cs(2)*exp(-(x0-cs(1)).^2/cs(2)^2/2)-y0;
cs0=rand(2,1);
cs=lsqnonlin(f,cs0);

所得结果为 $ \mu=0,\sigma=1$

lsqnonneg函数

求解非负的 $x$ ，使得
$$\underset{x}{\min }\Vert Cx-d{\Vert }_2^2$$
Matlab中的函数为
$$x=lsqnonneg(C,d,options)$$

例：已知 $C=\left[\begin{array}{cc}0.0372& 0.2869\\ 0.6861& 0.7071\\ 0.6233& 0.6245\\ 0.6344& 0.6170\end{array}\right]$ ，$d=\left[\begin{array}{c}0.8587\\ 0.1781\\ 0.0747\\ 0.8405\end{array}\right]$ ，求 $x=[x_1,x_2{]}^T$ 使得其满足 ${\min }_x\Vert Cx-d{\Vert }_2^2$ 。

c=[0.0372 0.2869;0.6861 0.7071;0.6233 0.6245;0.6344 0.6170];
d=[0.8587;0.1781;0.0747;0.8405];
x=lsqnonneg(c,d);

求得结果为 $x_1=0,x_2=0.6929$

Matlab中的曲线拟合用户图形界面解法

Matlab工具箱提供了命令cftool ，该命令给出了以为数据拟合的交互式环境。具体步骤如下

• 把数据导入到工作空间中
• 运行cftool ，打开用户图形界面窗口
• 选择适当的模型及逆行拟合
• 生成一些相关的统计量，并进行预测

曲线拟合与函数逼近

曲线拟合是已知一组离散数据 { ( x i , y i ) , i = 1 , ⋯   , n } \{(x_i,y_i),i=1,\cdots,n\} {(xi​,yi​),i=1,⋯,n} ，选择一个较简单的函数 $f(x)$ 在一定准则下，最接近这些数据。

如果已知一个较为复杂的连续函数 $y(x),x\in [a,b]$ ，要求选择一个较为简单的函数 $f(x)$ ，在一定准则下最接近 $y(x)$ ，就是所谓的函数逼近。

与曲线拟合的最小二乘准则相对应，函数逼近最常用的准则是最小平方逼近，即
$$J={\int }_a^b[f(x)-y(x){]}^{2}dx$$
达到最小，与函数拟合一样，选取一组函数 { r k ( x ) , k = 1 , ⋯   , m } \{r_k(x),k=1,\cdots,m\} {rk​(x),k=1,⋯,m} 构造函数 $f(x)$ ，即令
$$f(x)=a_1{r}_1(x)+\cdots +a_m{r}_m(x)$$
代入上式，求 $a_1,\cdots ,a_m$ 使得 $J$ 达到最小。利用极值必要条件可以得到

$$\left[\begin{array}{ccc}(r_1,r_1)& \cdots & (r_1,r_m)\\ ⋮& ⋮& ⋮\\ (r_m,r_1)& \cdots & (r_m,r_m)\end{array}\right]\left[\begin{array}{c}{a}_1\\ ⋮\\ a_m\end{array}\right]=\left[\begin{array}{c}(y,r_1)\\ ⋮\\ (y,r_m)\end{array}\right]$$

这里 $(g,h)={\int }_a^{b}g(x)h(x)dx$。当上述方程组系数矩阵行列式不为零的时候，有唯一解

最简单的是用多项式逼近函数，即选 $r_1(x)=1,r_2(x)=x,\cdots$。并且如果能使 ${\int }_a^b{r}_i(x)r_j(x)dx=0，i\ne j$，则方程组的系数矩阵将是对角阵，计算极大简化。满足这种性质的多项式称为正交多项式。

勒让德多项式就是在 $[-1,1]$ 区间上的正交多项式，它的表达式为
$$P_0(x)=1,P_k(x)=\frac{1}{2^{k}k!}\frac{d^k}{d{x}^k}(x^2-1{)}^k,k=1,2,\cdots 。$$
可以证明的是
$$\begin{gathered} {\int }_{-1}^1{P}_i(x)P_j(x)dx=\{\begin{array}{ll}0& i\ne j\\ \frac{2}{2i+1}& i=j\end{array} \\ {P}_{k+1}(x)=\frac{2k+1}{k+1}x{P}_k(x)-\frac{k}{k+1}{P}_{k-1}(x),k=1,2,\cdots \end{gathered}$$
常用的正交多项式还有第一类切比雪夫多项式
$$T_n(x)=\cos (narccosx),x\in [-1,1],n=0,1,2,\cdots$$
和拉盖尔多项式
$$L_n(x)=e^x\frac{d^n}{d{x}^n}(x^n{e}^{-x}),x\in [0,+\infty ],n=0,1,2,\cdots$$

例：求 $f(x)=cosx,x\in [-\frac{\pi }{2},\frac{\pi }{2}]$ 在 H = S p a n { 1 , x 2 , x 4 } H=Span\{1,x^2,x^4\} H=Span{1,x2,x4} 中的最佳平方逼近多项式

syms x
base=[1,x^2,x^4];
y1=base.' * base;
y2=cos(x)*base.';
r1=int(y1,-pi/2,pi/2);
r2=int(y2,-pi/2,pi/2);
a=r1\r2;
xs=double(a);

所得结果为 $xs=\left[\begin{array}{c}0.9996\\ -0.4964\\ 0.0372\end{array}\right]$ ，最佳平方逼近多项式为
$$y=0.9996-0.4964x^2+0.0372x^4$$

黄河小浪底调水调沙问题

问题的提出

2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄防水，形成人造洪峰进行调沙试验获得成功。整个试验期共20天，小浪底从6月19号开始预泄防水，直到7月13日结束并恢复正常供水。小浪底水利工程按设计拦沙量为75.5亿 $m^3$ 。在这之前，小浪底共积泥沙达 14,15亿吨。这次调水调沙试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙，在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于6月29日先后到达小浪底，7月3日达到最大流量2700$m^3/s$ ，使小浪底水库的排沙量不断增加。下表是试验数据。

现在根据试验数据，建立数学模型，研究以下问题：

(1) 给出估计任意时刻的排沙量以及总排沙量的方法

(2) 确定排沙量于水流量之间的关系

模型的建立与求解

已知给定的观测时间是等间隔的，每次的观测时刻可以记录为
$$t_i=3600(12i-4),i=1,\cdots ,24$$
第1次观测的时刻为 $t_1=28800$，最后一次观测的时刻为 $t_{24}=1022400$

记第 $i(i=1,2,\cdots ,24)$ 次观测时水流量为 $v_i$ ，含沙量为 $c_i$ ，则第 $i$ 次观测时的排沙量为 $y_i=c_i{v}_i$

对于问题 (1)，根据所给问题的试验数据，要计算任意时刻的排沙量，就要确定出排沙量随时间变化的规律，可以通过插值来实现。考虑到实际中的排沙量应该是时间的连续函数，为了提高模型的精度，可以采用三次样条函数进行插值。

利用Matlab函数，求出三次样条函数，得到排沙量 $y=y(x)$ 与时间的关系，然后进行积分，就可以得到总的排沙量
$$z={\int }_{t_1}^{t_{24}}y(t)dt$$
最后求得总的排沙量为 $1.844\times 10^{1}1t$ ，具体代码如下

%先把水流量和含沙量存到纯文本文件data.txt中
load data.txt
liu=data([1,3],:);
liu=liu';
liu=liu(:);
sha=data([2,4],:);
sha=sha';
sha=sha(:);
y=sha.*liu;
y=y';
i=1:24;
t=(12*i-4)*3600;
t1=t(1);te=t(end);
pp=csape(t,y);
ii=1:0.1:24;
ti=(12*ii-4)*3600;
t1=t(1);te=t(end);
pp=csape(t,y);
plot(t,y,'*');
hold on;
plot(ti,fnval(pp,ti));
TL=quadl(@(tt)fnval(pp,tt),t1,te);

对于问题 (2)，研究排沙量与水流量的关系，从试验数据可以看出，排沙量是随着水流量的增加而增长的，而后是随水流量的减少而减少的。显然，变化规律并非线性关系。因此，把问题分成两个部分，从开始水流量增加到最大值为第一阶段，从水流量的最导致到结束为第二阶段，分别来研究水流量和排沙量的关系。

绘制散点图

subplot(1,2,1),plot(liu(1:11),y(1:11),'*');
subplot(1,2,2),plot(liu(12:24),y(12:24),'*');

从散点图可以看出，第一阶段基本上符合线性关系。第一阶段和第二阶段都准备采用一次和二次曲线进行拟合，拿个模型的剩余标准差小，就采用哪个模型。

format long e
for j=1:2
    nihe1{j}=polyfit(liu(1:11),y(1:11),j);
    yhat1{j}=polyval(nihe1{j},liu(1:11));
    cha1(j)=sum((y(1:11)-yhat1{j}).^2);
    rmse1(j)=sqrt(cha1(j)/(10-j));
end
for j=1:2
    nihe2{j}=polyfit(liu(12:24),y(12:24),j+1);
    yhat2{j}=polyval(nihe2{j},liu(12:24));
    cha2(j)=sum((y(12:24)-yhat2{j}).^2);
    rmse2(j)=sqrt(cha1(j)/(11-j));
end
p1=nihe1{find(min(rmse1))};
p2=nihe2{find(min(rmse2))};
x1=min(liu(1:11)):max(liu(1:11));
x2=min(liu(12:24)):max(liu(12:24));
subplot(1,2,1),plot(liu(1:11),y(1:11),'*',x1,polyval(p1,x1));
subplot(1,2,2),plot(liu(12:24),y(12:24),'*',x2,polyval(p2,x2));

第一阶段 $y=250.5655v-373384.4661$

第二阶段 $y=0.1067v^2-180.4668+72421.0982$