插值与拟合
在实际问题中，一个函数
y
=
f
(
x
)
y=f(x)
y=f(x) 往往是通过实验观测得到的，仅已知函数
f
(
x
)
f(x)
f(x) 在某区间
[
a
,
b
]
[a,b]
[a,b] 撒谎给你一些列点的值，当需要这些节点
x
0
,
x
1
,
⋯
,
x
n
x_0,x_1,\cdots ,x_n
x0,x1,⋯,xn 之间的点上的函数值的时候，常用的简单的满足一定条件的函数
φ
(
x
)
\varphi (x)
φ(x) 替代
f
(
x
)
f(x)
f(x)，插值法就是一种常用的方法，其插值函数
φ
(
x
)
\varphi (x)
φ(x) 满足
φ
(
x
i
)
=
y
i
,
i
=
1
,
⋯
,
n
\varphi(x_i)=y_i, \quad i=1,\cdots,n
φ(xi)=yi,i=1,⋯,n
拟合也是已知有限个数据点，求近似函数，不要求过已知数据点，只要求在某种意义上他的这些点的总偏差最小
插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上是完全不同的。
插值方法
在平面上给定一组离散的点列，要求用一条曲线把这些点按照次序连接起来，称为插值
分段线性插值
简单来说，将没两个相邻的节点用直线连接起来，如此形成一条折线就是分段线性插值函数，记作 I n ( x ) I_n(x) In(x)，满足 I n ( x i ) = y i I_n(x_i)=y_i In(xi)=yi，且在每个小区间 [ x i , x i + 1 ] [x_i,x_{i+1}] [xi,xi+1] 上是线性函数
I
n
(
x
)
I_n(x)
In(x) 可以表示成
I
n
(
x
)
=
∑
i
=
0
n
y
i
l
i
(
x
)
I_n(x)=\sum\limits_{i=0}^{n}y_il_i(x)
In(x)=i=0∑nyili(x)，其中
l
i
=
{
x
−
x
i
−
1
x
i
−
x
i
−
1
x
∈
[
x
i
−
1
,
x
i
]
,
i
≠
0
x
−
x
i
+
1
x
i
−
x
i
+
1
x
∈
[
x
i
,
x
i
+
1
]
,
i
≠
n
0
l_i= \begin{cases} \frac{x-x_{i-1}}{x_i-x_{i-1}} \quad x \in [x_{i-1},x_i],i \neq 0 \\ \frac{x-x_{i+1}}{x_i-x_{i+1}} \quad x \in [x_i,x_{i+1}],i \neq n \\ 0 \end{cases}
li=⎩⎪⎨⎪⎧xi−xi−1x−xi−1x∈[xi−1,xi],i=0xi−xi+1x−xi+1x∈[xi,xi+1],i=n0
要求
I
n
(
x
)
I_n(x)
In(x) 具有良好的收敛性，即对于
x
∈
[
a
,
b
]
x \in [a,b]
x∈[a,b] 有
lim
x
→
∞
I
n
(
x
)
=
f
(
x
)
\lim _{x \to \infty} I_n(x)=f(x)
x→∞limIn(x)=f(x)
用
I
n
(
x
)
I_n(x)
In(x) 计算
x
x
x 点的插值时，只用到
x
x
x 左右两个节点，计算量与节点数
n
n
n 无关，但是
n
n
n 越大，分段越多，插值误差越小。实际上用函数表作插值计算时，分段线性插值就足够了。
拉格朗日插值多项式
拉格朗日插值多项式的基函数为
KaTeX parse error: \cr valid only within a tabular/array environment
l
i
(
x
)
l_i(x)
li(x) 是一个
n
n
n 次多项式，且满足
l
i
(
x
j
)
=
{
0
,
j
≠
i
1
,
j
=
i
l_i(x_j)= \begin{cases} 0,&j \neq i \\ 1,&j = i \end{cases}
li(xj)={0,1,j=ij=i
拉格朗日插值函数为
L
n
(
x
)
=
∑
i
=
0
n
y
i
l
i
(
x
)
=
∑
i
=
0
n
y
i
(
∏
j
=
0
j
≠
i
n
x
−
x
j
x
i
−
x
j
)
L_n(x)=\sum _{i=0} ^{n}y_il_i(x)=\sum _{i=0} ^{n}y_i(\prod _{j=0\\j \neq i} ^{n} \frac{x-x_j}{x_i-x_j})
Ln(x)=i=0∑nyili(x)=i=0∑nyi(j=0j=i∏nxi−xjx−xj)
样条插值
在许多工程技术中，提出的计算问题对插值函数的光滑性有较高的要求，要求曲线不仅要连续，还要有连续的曲率，这就产生了样条插值
-
样条插值的概念
样条本来是工程设计中的使用的一种绘图工具，是富有弹性的细木条或细金属条。绘图员利用它把一些已知的点连接成一条光滑的曲线，并使连接点处具有连续的曲率。三次样条插值就是由此抽象出来的。
数学上，将具有一定光滑性的分段多项式称为样条函数。具体说，给定区间 [ a , b ] [a,b] [a,b] 的一个分划
Δ : a = x 0 < x 1 < ⋯ < x n − 1 < x n = b \Delta:a=x_0<x_1<\cdots <x_{n-1}<x_n=b Δ:a=x0<x1<⋯<xn−1<xn=b
如果函数 S ( x ) S(x) S(x) 满足- 在每个小区间 [ x i , x i + 1 ] [x_i,x_{i+1}] [xi,xi+1] 上都是 S ( x ) S(x) S(x) 的 m m m 次多项式
- S ( x ) S(x) S(x) 在 [ a , b ] [a,b] [a,b] 上具有 m − 1 m-1 m−1 阶连续导数
则称 S ( x ) S(x) S(x) 为关于分划 Δ \Delta Δ 的 m m m 次样条函数，其图形为 m m m 次样条曲线
-
三次样条插值
利用样条函数进行插值，即插值函数为样条函数，称为样条插值。
下面介绍三次样条插值，即已知函数 y = f ( x ) y=f(x) y=f(x) 在区间 [ a , b ] [a,b] [a,b] 上的 n + 1 n+1 n+1 个节点
a = x 0 < x i < ⋯ < x n − 1 < x n = b a=x_0<x_i<\cdots<x_{n-1}<x_n=b a=x0<xi<⋯<xn−1<xn=b
上的值 y i = f ( x i ) y_i=f(x_i) yi=f(xi) ，求插值函数 S ( x ) S(x) S(x) ，使得- S ( x i ) = y i S(x_i)=y_i S(xi)=yi
- 在每个小区间 [ x i , x i + 1 ] [x_i,x_{i+1}] [xi,xi+1] 上， S ( x ) S(x) S(x) 都是三次多项式，记为 S i ( x ) S_i(x) Si(x)
- S ( x ) S(x) S(x) 在 [ a , b ] [a,b] [a,b] 上二阶连续可微
这样的函数 S ( x ) S(x) S(x) 称为 f ( x ) f(x) f(x) 的三次样条插值函数
由条件2，不妨记
S ( x ) = S i ( x ) , x ∈ [ x i , x i + 1 ] , i = 0 , 1 , ⋯ , n − 1 S i ( x ) = a i x 3 + b i x 2 + c i x + d \begin{aligned} &S(x)={S_i(x),x \in [x_i,x_{i+1}],i=0,1,\cdots ,n-1} \\ &S_i(x)=a_ix^3+b_ix^2+c_ix+d \end{aligned} S(x)=Si(x),x∈[xi,xi+1],i=0,1,⋯,n−1Si(x)=aix3+bix2+cix+d
式*有 4 n 4n 4n 个待定参数由条件3，不妨记
{ S i ( x i + 1 ) = S i + 1 ( x i + 1 ) S i ′ ( x i + 1 ) = S i + 1 ′ ( x i + 1 ) S i ′ ′ ( x i + 1 ) = S i + 1 ′ ′ ( x i + 1 ) \begin{cases} S_i(x_{i+1})=S_{i+1}(x_{i+1}) \\ S_i'(x_{i+1})=S_{i+1}'(x_{i+1}) \\ S_i''(x_{i+1})=S_{i+1}''(x_{i+1}) \end{cases} ⎩⎪⎨⎪⎧Si(xi+1)=Si+1(xi+1)Si′(xi+1)=Si+1′(xi+1)Si′′(xi+1)=Si+1′′(xi+1)
结合条件1和条件3，共能够得到 4 n − 2 4n-2 4n−2 个方程，为了求出条件2的待定参数，还需要由两个边界条件常用的三次样条函数的边界条件由3中类型
-
S ′ ( a ) = y 0 ′ , S ′ ( b ) = y n ′ S'(a)=y'_0,S'(b)=y'_n S′(a)=y0′,S′(b)=yn′。
由这种边界条件建立的样条插值函数称为 f ( x ) f(x) f(x) 的完备三次样条插值函数。
特别地，当 y 0 ′ = y n ′ = 0 y'_0=y'_n=0 y0′=yn′=0 时，样条曲线在端点处呈现水平状态
如果 f ′ ( x ) f'(x) f′(x) 不知道，我们可以要求 S ′ ( x ) S'(x) S′(x) 与 f ′ ( x ) f'(x) f′(x) 在端点处近似相等，这时以 x 0 , x 1 , x 2 , x 3 x_0,x_1,x_2,x_3 x0,x1,x2,x3 为节点做一个三次 Newton 插值多项式 N a ( x ) N_a(x) Na(x)，以 x n , x n − 1 , x n − 2 , x n − 3 x_n,x_{n-1},x_{n-2},x_{n-3} xn,xn−1,xn−2,xn−3 做一个三次 Newton 插值多项式 N b ( x ) N_b(x) Nb(x) ，要求
S ′ ( a ) = N ′ ( a ) , S ′ ( b ) = N ′ ( b ) S'(a)=N'(a),S'(b)=N'(b) S′(a)=N′(a),S′(b)=N′(b)
由这种边界条件建立的三次样条称为 f ( x ) f(x) f(x) 的Lagrange三次样条插值函数 -
S ′ ′ ( a ) = y 0 ′ ′ , S ′ ′ ( b ) = y n ′ ′ S''(a)=y''_0,S''(b)=y''_n S′′(a)=y0′′,S′′(b)=yn′′
特别地，当 y 0 ′ ′ = y n ′ ′ y''_0=y''_n y0′′=yn′′ 时，称为自然边界条件
-
S ′ ( a + 0 ) = S ′ ( b − 0 ) , S ′ ′ ( a + 0 ) = S ′ ′ ( b − 0 ) S'(a+0)=S'(b-0),S''(a+0)=S''(b-0) S′(a+0)=S′(b−0),S′′(a+0)=S′′(b−0)
此条件称为周期条件
Matlab插值工具箱
-
一维插值函数
Matlab中由现成的一维插值函数
interp1，其基本语法为
y = i n t e r p 1 ( x 0 , y 0 , x , ′ m e t h o d ′ ) y=interp1(x_0,y_0,x,'method') y=interp1(x0,y0,x,′method′)
其中method指定插值的方法，其默认为线性插值。其全部参数如下参数 方法 nearest最近项插值 linear线性插值 spline立方样条 cubic立方插值 所有的插值方法都要求 x 0 x_0 x0 是单调的
当 x 0 x_0 x0 为等距时可以用快速插值法，使用快速插值法的格式为
*nearest``*linear``*spline``*cubic -
三次样条插值
在Matlab中数据点称为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非纽结条件。这种条件强迫第1个和第2个三次多项式的三阶导数相等。对最后一个和倒数第二个两个三次多项式也做同样的处理。
在Matlab中的三次样条插值有一下函数
y = i n t e r p 1 ( x 0 , y 0 , x , ′ s p l i n e ′ ) y = s p l i n e ( x 0 , y 0 , x ) p p = c s a p e ( x 0 , y 0 , x ) o r p p = c s a p e ( x 0 , y 0 , c o n d s , v a l c o n d s ) y = f n v a l ( p p , x ) \begin{aligned} &y=interp1(x_0,y_0,x,'spline') \\ &y=spline(x_0,y_0,x) \\ &pp=csape(x_0,y_0,x) \ \ or \ \ pp=csape(x_0,y_0,conds,valconds) \\ &y=fnval(pp,x) \end{aligned} y=interp1(x0,y0,x,′spline′)y=spline(x0,y0,x)pp=csape(x0,y0,x) or pp=csape(x0,y0,conds,valconds)y=fnval(pp,x)
其中： x 0 , y 0 x_0,y_0 x0,y0 为已知的数据点， x , y x,y 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 ] [0.0] [0.0] -
variational设置边界的二阶导数值为 [ 0 , 0 ] [0,0] [0,0]
对于一些特殊的边界条件，可以通过
conds的一个 1 × 2 1 \times 2 1×2 矩阵来表示，conds元素的取值为0，1，2。conds(i)=j的含义是给定端点 i i i 的 j j j 阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件，例：conds=[2,1]表示左边界为2阶导数，右边界为1阶导数，其值由valconds参数给出例：待加工零件的外形根据工艺要求由一组数据 ( x , y ) (x,y) (x,y) 给出，加工时每一刀只能沿 x x x 轴或者 y y y 轴方向走非常小一步，这就需要从已知的数据得到加工所要求的步长的很小的 ( x , y ) (x,y) (x,y) 坐标
表中给出的 x , y x,y x,y 数据位于机翼断面的下轮廓线上，假设需要得到 x x x 坐标的每改变0.1时的 y y y 坐标。试完成加工所需要的数据，画出曲线，并求出 x = 0 x=0 x=0 处的曲线斜率和 13 ≤ x ≤ 15 13 \leq x \leq 15 13≤x≤15 范围内 y y y 的最小值。要求分别采用分段线性和三次样条两种插值方法
x x x 0 3 5 7 9 11 12 13 14 15 y y 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]结果如图
-
-
二维插值
前面讲述的都是一维插值，即节点为一维变量，插值函数时一元函数。如果节点为二维，插值函数就是二元函数。如果在某区域测量若干节点的高程，为了画出较为精确的等高线图，要先插入更多的点，计算这些点的高程，这就是二维插值。
1）插值节点为网格节点
已知 m × n m \times n m×n 个节点： ( x i , y j , z i j ) (x_i,y_j,z_{ij}) (xi,yj,zij) ，且满足 x 1 < ⋯ < x m ; y i < ⋯ < y n x_1 < \cdots <x_m;y_i < \cdots < y_n x1<⋯<xm;yi<⋯<yn 。求在 ( x , y ) (x,y) (x,y) 处的插值 z z z。
Matlab中有一些计算二维插值的命令
z = i n t e r p 2 ( x 0 , y 0 , z 0 , x , y , ′ m e t h o d ′ ) z=interp2(x0,y0,z0,x,y,'method') z=interp2(x0,y0,z0,x,y,′method′)
其中： x 0 , y 0 , z 0 x0,y0,z0 x0,y0,z0 分别为 m m m 维和 n n n 维向量，表示节点； z 0 z0 z0 为 n × m n\times m n×m 矩阵，表示节点值； x , y x,y x,y 为一维数组，表示插值点， x x x 与 y y y 应是方向不同的向量，即一个是行向量一个是列向量。 z z z 是矩阵，行数为 y y y 的维数，列数为 x x 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 x x 和 y y y 方向每隔100 m m 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 ) (166,177) (166,177)，高程为 z = 720.7540 z=720.7540 z=720.7540。
2）插值节点为散乱节点
已知 n n n 个节点 ( x i , y i , z i ) (x_i,y_i,z_i) (xi,yi,zi) ，求点 ( x , y ) (x,y) (x,y) 处的插值 z z z 。
对于上述问题，Matlab中提供了插值函数
griddata
Z I = g r i d d a t a ( x , y , z , X I , Y I , ′ m e t h o d ′ ) ZI=griddata(x,y,z,XI,YI,'method') ZI=griddata(x,y,z,XI,YI,′method′)
其中， x , y , z x,y,z x,y,z 均为 n n n 维向量，指出所给数据点的横坐标、纵坐标和竖坐标。向量 X I , Y I XI,YI XI,YI 是给定网格点的横坐标和纵坐标。返回值 Z I ZI ZI 为网格 ( X I , Y I ) (XI,YI) (XI,YI) 处的函数值。 X I XI XI 为行向量， Y I YI YI 为列向量参数 方法 cubic立方插值 nearest最近点插值 P.S. Matlab插值时外插值有时是不确定的，可以采用立方插值和最近点插值混合插值法，把不确定的插值换成最近点插值。
例：在某海域内测量得到一些点 ( x , y ) (x,y) (x,y) 处的水深 z z 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 n n 个点 ( x i , y i ) (x_i,y_i) (xi,yi) ， x i x_i xi 互不相同，求一个函数 y = f ( x ) y=f(x) y=f(x) ，使 f ( x ) f(x) f(x) 在某种准则下与所有数据点最为接近，即曲线拟合效果最好
线性最小二乘法使解决曲线拟合的最常用的方法，其基本思路为：令
f
(
x
)
=
a
1
r
1
(
x
)
+
a
2
r
2
(
x
)
+
⋯
+
a
m
r
m
(
x
)
f(x)=a_1r_1(x)+a_2r_2(x)+\cdots +a_mr_m(x)
f(x)=a1r1(x)+a2r2(x)+⋯+amrm(x)
式中：
r
k
(
x
)
r_k(x)
rk(x) 为事先选定的一组线性无关的函数；
a
k
a_k
ak 为待定系数
拟合准则是使 y i ( i = 1 , 2 , ⋯ , n ) y_i(i=1,2,\cdots ,n) yi(i=1,2,⋯,n) 与 f ( x i ) f(x_i) f(xi) 的距离 δ i \delta_i δi 的平方和最小，称为最小二乘准则。
-
系数 a k a_k ak 的确定
记
J ( a 1 , ⋯ , a m ) = ∑ i = 1 n δ i 2 = ∑ i = 1 n [ f ( x i ) − y i ] 2 J(a_1,\cdots ,a_m)=\sum _{i=1} ^{n} \delta _i ^2=\sum _{i=1} ^{n} [f(x_i)-y_i]^2 J(a1,⋯,am)=i=1∑nδi2=i=1∑n[f(xi)−yi]2
为求 a 1 , ⋯ , a m a_1,\cdots,a_m a1,⋯,am 使 J J J 达到最小，只需要利用极值的必要条件 KaTeX parse error: Undefined control sequence: \part at position 7: \frac{\̲p̲a̲r̲t̲ ̲J}{\part a_j}=0… ，得到关于 a 1 , ⋯ , a m a_1,\cdots ,a_m a1,⋯,am 的线性方程组
∑ i = 1 n r j ( x i ) [ ∑ k = 1 m a k r k ( x i ) − y i ] = 0 , j = 1 , ⋯ , m \sum _{i=1} ^{n} r_j(x_i)[\sum _{k=1} ^ma_kr_k(x_i)-y_i]=0,\quad j=1,\cdots ,m i=1∑nrj(xi)[k=1∑makrk(xi)−yi]=0,j=1,⋯,m
即
∑ k = 1 m a k [ ∑ i = 1 n r j ( x i ) r k ( x i ) ] = ∑ i = 1 n r j ( x i ) y i j = 1 , ⋯ , 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 \quad j=1,\cdots,m k=1∑mak[i=1∑nrj(xi)rk(xi)]=i=1∑nrj(xi)yij=1,⋯,m
记
R = [ r 1 ( x 1 ) ⋯ r m ( x 1 ) ⋮ ⋮ r 1 ( x n ) ⋯ r m ( x n ) ] A = [ a i , ⋯ , a m ] T Y = [ y 1 , ⋯ , y n ] T \mathbf{R}= \left[\begin{matrix} r_1(x_1)&\cdots&r_m(x_1)\\ \vdots &&\vdots \\ r_1(x_n)&\cdots&r_m(x_n) \end{matrix}\right] \\ \mathbf{A}=[a_i,\cdots,a_m]^T\\ \mathbf{Y}=[y_1,\cdots,y_n]^T R=⎣⎢⎡r1(x1)⋮r1(xn)⋯⋯rm(x1)⋮rm(xn)⎦⎥⎤A=[ai,⋯,am]TY=[y1,⋯,yn]T
由此上述方程组可以写成
R T R A = R T Y \mathbf{R}^T \mathbf{R} \mathbf{A}=\mathbf{R}^T\mathbf{Y} RTRA=RTY
当 { r 1 ( x ) ⋯ , r m ( x ) } \{r_1(x)\,\cdots,r_m(x)\} {r1(x)⋯,rm(x)} 线性无关的时候， R \mathbf{R} R 列满秩， R T R \mathbf{R^TR} RTR 可逆，于是方程组有唯一解
A = ( R T R ) − 1 R T Y \mathbf{A}=\mathbf{(R^TR)^{-1}R^TY} A=(RTR)−1RTY -
函数 r k ( x ) r_k(x) rk(x) 的选取
面对一组数据 ( x i , y i ) (x_i,y_i) (xi,yi) 用线性最小二乘法做曲线拟合的时候，最关键的一步就是恰当地选取 r 1 ( x ) , ⋯ , r m ( x ) r_1(x),\cdots,r_m(x) r1(x),⋯,rm(x) 。如果通过机理分析，能够知道 y y y 与 x x x 之间的函数关系，则 r 1 ( x ) , ⋯ , r m ( x ) r_1(x),\cdots,r_m(x) r1(x),⋯,rm(x) 容易确定。若无法知道 y y y 与 x x x 之间的函数关系，通常可以通过将数据 ( x i , y i ) (x_i,y_i) (xi,yi) 作图，直观判断应该用什么曲线去拟合。通常有一下形式
- 直线 y = a x + b y=ax+b y=ax+b
- 多项式 y = a m x m + ⋯ + a 1 x + a 0 y=a_mx^m+\cdots+a_1x+a_0 y=amxm+⋯+a1x+a0
- 双曲线 y = a 1 x + a 2 y=\frac{a_1}{x}+a_2 y=xa1+a2
- 指数曲线 y = a 1 e a 2 x y=a_1e^{a_2x} y=a1ea2x （对于指数曲线，拟合前要取对数进行变量代换转化成线性函数）
最小二乘法的Matlab实现
-
解方程组方法
在上面的记号下，
J ( a 1 , ⋯ , a m ) = ∥ R A − Y ∥ 2 2 J(a_1,\cdots,a_m)=\|\mathbf{RA-Y}\|_2^2 J(a1,⋯,am)=∥RA−Y∥22
Matlab中的线性最小二乘的标准型是
min A ∥ R A − Y ∥ 2 2 \min _A \|\mathbf{RA-Y}\|_2^2 Amin∥RA−Y∥22
系数矩阵 A = R \ Y \mathbf{A}=\mathbf{R} \backslash \mathbf{Y} A=R\Y其中： R \mathbf{R} R 为方程组的系数矩阵
例：用最小二乘法求一个形如 y = a + b x 2 y=a+bx^2 y=a+bx2 的经验公式，使其与下表拟合
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.05 x 2 y=0.9726+0.05x^2 y=0.9726+0.05x2
-
多项式拟合方法
如果选取 { 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 m m 次多项式拟合给定数据，Matlab中现有的函数为
polyfit，具体命令如下
p = p o l y f i t ( x 0 , y 0 , m ) p=polyfit(x0,y0,m) p=polyfit(x0,y0,m)
其中： x 0 , y 0 x0,y0 x0,y0 为要拟合的数据， m m m 为多项式的次数，返回值 p p p 为多项式的系数向量多项式在 x x x 处的值可以用函数
polyval来计算，具体命令如下
y = p o l y v a l ( p , x ) y=polyval(p,x) 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 y_{2017}=233.4286,y_{2018}=253.9286 y2017=233.4286,y2018=253.9286
最小二乘优化
在无约束最优化问题中，有些重要的特殊情形，比如目标函数是由若干个函数的平方和构成的，这类函数一般可以写成
F
(
x
)
=
∑
i
=
1
m
f
i
2
(
x
)
,
x
∈
R
n
F(x)=\sum _{i=1} ^{m} f_i^2(x),\quad x \in R^n
F(x)=i=1∑mfi2(x),x∈Rn
式中：
x
=
[
x
1
,
⋯
,
x
n
]
T
x=[x_1,\cdots ,x_n]^T
x=[x1,⋯,xn]T
把极小化这类函数的问题
min
F
(
x
)
=
∑
i
=
1
m
f
i
2
(
x
)
\min F(x)=\sum _{i=1} ^{m} f_i^2(x)
minF(x)=i=1∑mfi2(x)
称为最小二乘优化问题
最小二乘优化问题是一类比较特殊的优化问题，在处理这类问题的时候，Matlab也提供了一些强大的函数。在Matlab优化工具箱中，用于求解最小二乘优化问题的函数有lsqlin,lsqcurvefit,lsqnonlin,lsqnonneg
lsqlin函数
求解
min
x
1
2
∥
C
⋅
x
−
d
∥
2
2
s
.
t
.
{
A
⋅
x
≤
b
A
e
q
⋅
x
=
b
e
q
l
b
≤
x
≤
u
b
\min _x \frac{1}{2}\|C \cdot x-d\|_2 ^2 \\ s.t. \begin{cases} A \cdot x \leq b \\ A_{eq} \cdot x = b_{eq} \\ lb \leq x \leq ub \end{cases}
xmin21∥C⋅x−d∥22s.t.⎩⎪⎨⎪⎧A⋅x≤bAeq⋅x=beqlb≤x≤ub
式中：
C
,
A
,
A
e
q
C,A,A_{eq}
C,A,Aeq 均为矩阵，
d
,
b
,
b
e
q
,
l
b
,
u
b
d,b,b_{eq},lb,ub
d,b,beq,lb,ub 为向量
在Matlab中的命令为
x
=
l
s
q
l
i
n
(
C
,
d
,
A
,
b
,
A
e
q
,
b
e
q
,
l
b
,
u
b
,
x
0
)
x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)
x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)
例：用最小二乘法求一个形如 y = a + b x 2 y=a+bx^2 y=a+bx2 的经验公式，使其与下表拟合
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.05 x 2 y=0.9726+0.05x^2 y=0.9726+0.05x2
lsqcurvefit函数
给定输入输出数列
x
d
a
t
a
,
y
d
a
t
a
xdata,ydata
xdata,ydata 求参量
x
x
x ，使得
min
x
∥
F
(
x
,
x
d
a
t
a
)
−
y
d
a
t
a
∥
2
2
=
∑
i
(
F
(
x
,
x
d
a
t
a
i
)
−
y
d
a
t
a
i
)
2
\min _x \|F(x,xdata)-ydata\| _2 ^2=\sum _i (F(x,xdata_i)-ydata_i)^2
xmin∥F(x,xdata)−ydata∥22=i∑(F(x,xdatai)−ydatai)2
Matlab中的函数为
x
=
l
s
q
c
u
r
v
e
f
i
t
(
f
u
n
,
x
0
,
x
d
a
t
a
,
y
d
a
t
a
,
l
b
,
u
b
,
o
p
t
i
o
n
s
)
x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)
x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)
其中：
f
u
n
fun
fun 为定义函数
F
(
x
,
x
d
a
t
a
)
F(x,xdata)
F(x,xdata) 的.m文件
例：根据下列数据，拟合函数 y = e − k 1 x 1 s i n ( k 2 x 2 ) + x 3 2 y=e^{-k_1x_1}sin(k_2x_2)+x_3^2 y=e−k1x1sin(k2x2)+x32 中的参数 k 1 , k 2 k_1,k_2 k1,k2
y / k g y/kg y/kg x 1 / c m 2 x_1/cm^2 x1/cm2 x 2 / k g x_2/kg x2/kg x 3 / k g x_3/kg x3/kg y / k g y/kg y/kg x 1 / c m 2 x_1/cm^2 x1/cm2 x 2 / k g x_2/kg x2/kg x 3 / k g x_3/kg x3/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
)
,
⋯
,
f
k
(
x
)
]
T
F(x)=[f_1(x),\cdots,f_k(x)]^T
F(x)=[f1(x),⋯,fk(x)]T ，求
x
x
x 使得
min
x
∥
F
(
x
)
∥
2
2
\min _x \|F(x)\|_2 ^2
xmin∥F(x)∥22
在Matlab中的函数为
x
=
l
s
q
n
o
n
l
i
n
(
f
u
n
,
x
0
,
l
b
,
u
b
,
o
p
t
i
o
n
s
)
x=lsqnonlin(fun,x0,lb,ub,options)
x=lsqnonlin(fun,x0,lb,ub,options)
其中
f
u
n
fun
fun 为定义向量函数
F
(
x
)
F(x)
F(x) 的 .m 文件
例：用
lsqnonlin函数拟合 y = 1 2 π σ e − ( x − μ ) 2 2 σ 2 y=\frac{1}{\sqrt{2\pi}\sigma}e^{-\frac{(x- \mu)^2}{2\sigma^2}} y=2π σ1e−2σ2(x−μ)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
x
x ，使得
min
x
∥
C
x
−
d
∥
2
2
\min _x \|Cx-d\|_2^2
xmin∥Cx−d∥22
Matlab中的函数为
x
=
l
s
q
n
o
n
n
e
g
(
C
,
d
,
o
p
t
i
o
n
s
)
x=lsqnonneg(C,d,options)
x=lsqnonneg(C,d,options)
例：已知 C = [ 0.0372 0.2869 0.6861 0.7071 0.6233 0.6245 0.6344 0.6170 ] C=\left[ \begin{matrix} 0.0372&0.2869\\0.6861&0.7071\\0.6233&0.6245\\0.6344&0.6170 \end{matrix}\right] C=⎣⎢⎢⎡0.03720.68610.62330.63440.28690.70710.62450.6170⎦⎥⎥⎤ ， d = [ 0.8587 0.1781 0.0747 0.8405 ] d=\left[\begin{matrix}0.8587\\0.1781\\0.0747\\0.8405\end{matrix}\right] d=⎣⎢⎢⎡0.85870.17810.07470.8405⎦⎥⎥⎤ ，求 x = [ x 1 , x 2 ] T x=[x_1,x_2]^T x=[x1,x2]T 使得其满足 min x ∥ C x − d ∥ 2 2 \min _x \|Cx-d\|_2^2 minx∥Cx−d∥22 。
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 x_1=0,x_2=0.6929 x1=0,x2=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 ) f(x) f(x) 在一定准则下，最接近这些数据。
如果已知一个较为复杂的连续函数 y ( x ) , x ∈ [ a , b ] y(x),x\in [a,b] y(x),x∈[a,b] ，要求选择一个较为简单的函数 f ( x ) f(x) f(x) ，在一定准则下最接近 y ( x ) y(x) y(x) ，就是所谓的函数逼近。
与曲线拟合的最小二乘准则相对应，函数逼近最常用的准则是最小平方逼近，即
J
=
∫
a
b
[
f
(
x
)
−
y
(
x
)
]
2
d
x
J=\int_a ^b [f(x)-y(x)]^2dx
J=∫ab[f(x)−y(x)]2dx
达到最小，与函数拟合一样，选取一组函数
{
r
k
(
x
)
,
k
=
1
,
⋯
,
m
}
\{r_k(x),k=1,\cdots,m\}
{rk(x),k=1,⋯,m} 构造函数
f
(
x
)
f(x)
f(x) ，即令
f
(
x
)
=
a
1
r
1
(
x
)
+
⋯
+
a
m
r
m
(
x
)
f(x)=a_1r_1(x)+\cdots+a_mr_m(x)
f(x)=a1r1(x)+⋯+amrm(x)
代入上式，求
a
1
,
⋯
,
a
m
a_1,\cdots ,a_m
a1,⋯,am 使得
J
J
J 达到最小。利用极值必要条件可以得到
[ ( r 1 , r 1 ) ⋯ ( r 1 , r m ) ⋮ ⋮ ⋮ ( r m , r 1 ) ⋯ ( r m , r m ) ] [ a 1 ⋮ a m ] = [ ( y , r 1 ) ⋮ ( y , r m ) ] \left[ \begin{matrix} (r_1,r_1)&\cdots&(r_1,r_m)\\ \vdots&\vdots&\vdots\\ (r_m,r_1)&\cdots&(r_m,r_m) \end{matrix}\right] \left[\begin{matrix} a_1\\ \vdots\\ a_m \end{matrix}\right] = \left[\begin{matrix} (y,r_1)\\ \vdots\\ (y,r_m) \end{matrix}\right] ⎣⎢⎡(r1,r1)⋮(rm,r1)⋯⋮⋯(r1,rm)⋮(rm,rm)⎦⎥⎤⎣⎢⎡a1⋮am⎦⎥⎤=⎣⎢⎡(y,r1)⋮(y,rm)⎦⎥⎤
这里 ( g , h ) = ∫ a b g ( x ) h ( x ) d x (g,h)=\int _a ^b g(x)h(x)dx (g,h)=∫abg(x)h(x)dx。当上述方程组系数矩阵行列式不为零的时候，有唯一解
最简单的是用多项式逼近函数，即选 r 1 ( x ) = 1 , r 2 ( x ) = x , ⋯ r_1(x)=1,r_2(x)=x,\cdots r1(x)=1,r2(x)=x,⋯。并且如果能使 ∫ a b r i ( x ) r j ( x ) d x = 0 ， i ≠ j \int _a ^b r_i(x)r_j(x)dx=0，i\neq j ∫abri(x)rj(x)dx=0，i=j，则方程组的系数矩阵将是对角阵，计算极大简化。满足这种性质的多项式称为正交多项式。
勒让德多项式就是在
[
−
1
,
1
]
[-1,1]
[−1,1] 区间上的正交多项式，它的表达式为
P
0
(
x
)
=
1
,
P
k
(
x
)
=
1
2
k
k
!
d
k
d
x
k
(
x
2
−
1
)
k
,
k
=
1
,
2
,
⋯
。
P_0(x)=1,P_k(x)=\frac{1}{2^kk!}\frac{d^k}{dx^k}(x^2-1)^k,k=1,2,\cdots。
P0(x)=1,Pk(x)=2kk!1dxkdk(x2−1)k,k=1,2,⋯。
可以证明的是
∫
−
1
1
P
i
(
x
)
P
j
(
x
)
d
x
=
{
0
i
≠
j
2
2
i
+
1
i
=
j
P
k
+
1
(
x
)
=
2
k
+
1
k
+
1
x
P
k
(
x
)
−
k
k
+
1
P
k
−
1
(
x
)
,
k
=
1
,
2
,
⋯
\int _{-1} ^{1} P_i(x)P_j(x)dx= \begin{cases} 0&i\neq j\\ \frac{2}{2i+1}&i=j \end{cases}\\ P_{k+1}(x)=\frac{2k+1}{k+1}xP_k(x)-\frac{k}{k+1}P_{k-1}(x),k=1,2,\cdots
∫−11Pi(x)Pj(x)dx={02i+12i=ji=jPk+1(x)=k+12k+1xPk(x)−k+1kPk−1(x),k=1,2,⋯
常用的正交多项式还有第一类切比雪夫多项式
T
n
(
x
)
=
cos
(
n
a
r
c
c
o
s
x
)
,
x
∈
[
−
1
,
1
]
,
n
=
0
,
1
,
2
,
⋯
T_n(x)=\cos(narccosx),x \in [-1,1] ,n=0,1,2,\cdots
Tn(x)=cos(narccosx),x∈[−1,1],n=0,1,2,⋯
和拉盖尔多项式
L
n
(
x
)
=
e
x
d
n
d
x
n
(
x
n
e
−
x
)
,
x
∈
[
0
,
+
∞
]
,
n
=
0
,
1
,
2
,
⋯
L_n(x)=e^x\frac{d^n}{dx^n}(x^ne^{-x}),x \in [0,+\infin],n=0,1,2,\cdots
Ln(x)=exdxndn(xne−x),x∈[0,+∞],n=0,1,2,⋯
例：求 f ( x ) = c o s x , x ∈ [ − π 2 , π 2 ] f(x)=cosx,x \in [-\frac{\pi}{2},\frac{\pi}{2}] f(x)=cosx,x∈[−2π,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);所得结果为 x s = [ 0.9996 − 0.4964 0.0372 ] xs=\left[\begin{matrix}0.9996\\-0.4964\\0.0372\end{matrix}\right] xs=⎣⎡0.9996−0.49640.0372⎦⎤ ，最佳平方逼近多项式为
y = 0.9996 − 0.4964 x 2 + 0.0372 x 4 y=0.9996-0.4964x^2+0.0372x^4 y=0.9996−0.4964x2+0.0372x4
黄河小浪底调水调沙问题
问题的提出
2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄防水，形成人造洪峰进行调沙试验获得成功。整个试验期共20天，小浪底从6月19号开始预泄防水，直到7月13日结束并恢复正常供水。小浪底水利工程按设计拦沙量为75.5亿 m 3 m^3 m3 。在这之前，小浪底共积泥沙达 14,15亿吨。这次调水调沙试验一个重要目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙，在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于6月29日先后到达小浪底，7月3日达到最大流量2700 m 3 / s m^3/s m3/s ，使小浪底水库的排沙量不断增加。下表是试验数据。
现在根据试验数据，建立数学模型，研究以下问题：
(1) 给出估计任意时刻的排沙量以及总排沙量的方法
(2) 确定排沙量于水流量之间的关系
模型的建立与求解
已知给定的观测时间是等间隔的，每次的观测时刻可以记录为
t
i
=
3600
(
12
i
−
4
)
,
i
=
1
,
⋯
,
24
t_i=3600(12i-4),\quad i=1,\cdots,24
ti=3600(12i−4),i=1,⋯,24
第1次观测的时刻为
t
1
=
28800
t_1=28800
t1=28800，最后一次观测的时刻为
t
24
=
1022400
t_{24}=1022400
t24=1022400
记第 i ( i = 1 , 2 , ⋯ , 24 ) i(i=1,2,\cdots,24) i(i=1,2,⋯,24) 次观测时水流量为 v i v_i vi ，含沙量为 c i c_i ci ，则第 i i i 次观测时的排沙量为 y i = c i v i y_i=c_iv_i yi=civi
对于问题 (1)，根据所给问题的试验数据，要计算任意时刻的排沙量，就要确定出排沙量随时间变化的规律，可以通过插值来实现。考虑到实际中的排沙量应该是时间的连续函数，为了提高模型的精度，可以采用三次样条函数进行插值。
利用Matlab函数，求出三次样条函数，得到排沙量
y
=
y
(
x
)
y=y(x)
y=y(x) 与时间的关系，然后进行积分，就可以得到总的排沙量
z
=
∫
t
1
t
24
y
(
t
)
d
t
z=\int _{t_1} ^{t_{24}} y(t)dt
z=∫t1t24y(t)dt
最后求得总的排沙量为
1.844
×
1
0
1
1
t
1.844 \times10^11t
1.844×1011t ，具体代码如下
%先把水流量和含沙量存到纯文本文件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.5655 v − 373384.4661 y=250.5655v-373384.4661 y=250.5655v−373384.4661
第二阶段 y = 0.1067 v 2 − 180.4668 + 72421.0982 y=0.1067v^2-180.4668+72421.0982 y=0.1067v2−180.4668+72421.0982