插值
多项式插值
给定 n + 1 n+1 n+1 个互异点 a ⩽ x 0 < x 1 . . . < x n ⩽ b a \leqslant {x_0} < {x_1}... < {x_n} \leqslant b a⩽x0<x1...<xn⩽b 的值 y 0 , y 1 , ⋯ , y n y_0, y_1 , \cdots , y_n y0,y1,⋯,yn, 求一个次数不超过 n n n 次的多项式 p ( x ) = a 0 + a 1 x + . . . + a n x n p(x) = {a_0} + {a_1}x + ... + {a_n}{x^n} p(x)=a0+a1x+...+anxn, s . t . s.t. s.t.
p ( x i ) = y i ( i = 0 , … n ) p(x_i) = y_i (i = 0, … n) p(xi)=yi(i=0,…n)
称 p ( x ) p(x) p(x) 为 f ( x ) f(x) f(x) 的 n n n 次插值多项式
Lagrange 插值
l k ( x ) = ( x − x 0 ) ( x − x 1 ) ⋯ ( x − x k − 1 ) ( x − x k + 1 ) ⋯ ( x − x n ) ( x k − x 0 ) ( x k − x 1 ) ⋯ ( x k − x k − 1 ) ( x k − x k + 1 ) ⋯ ( x k − x n ) l_{k}(x)=\frac{\left(x-x_{0}\right)\left(x-x_{1}\right) \cdots\left(x-x_{k-1}\right)\left(x-x_{k+1}\right) \cdots\left(x-x_{n}\right)}{\left(x_{k}-x_{0}\right)\left(x_{k}-x_{1}\right) \cdots\left(x_{k}-x_{k-1}\right)\left(x_{k}-x_{k+1}\right) \cdots\left(x_{k}-x_{n}\right)} lk(x)=(xk−x0)(xk−x1)⋯(xk−xk−1)(xk−xk+1)⋯(xk−xn)(x−x0)(x−x1)⋯(x−xk−1)(x−xk+1)⋯(x−xn)
即
∏ i ≠ k x − x i x k − x i \prod_{{ i \neq k}} \frac{x-x_{i}}{x_{k}-x_{i}} i=k∏xk−xix−xi
k = 0 , 1 , ⋯ , n k=0, 1, \cdots , n k=0,1,⋯,n
满足
l k ( x i ) = δ k i = { 1 , i = k 0 , i ≠ k {l_k}({x_i}) = {\delta _{ki}} = \left\{ \begin{gathered} 1,{\text{ }}i = k \\ 0,{\text{ }}i \ne k \\ \end{gathered} \right. lk(xi)=δki={1, i=k0, i=k
l 0 ( x ) , l 1 ( x ) , ⋯ , l n ( x ) l_0(x), l_1(x), \cdots , l_n(x) l0(x),l1(x),⋯,ln(x) 线性无关,即它们构成线性空间 P n ( x ) P_n(x) Pn(x) 的一组基
则
p ( x ) = y 0 l 0 ( x ) + y 1 l 1 ( x ) + ⋯ + y n l n ( x ) p(x)=y_{0} l_{0}(x)+y_{1} l_{1}(x)+\cdots+y_{n} l_{n}(x) p(x)=y0l0(x)+y1l1(x)+⋯+ynln(x)
满足插值条件 p ( x i ) = y i ( i = 0 , … n ) p(x_i) = y_i \ (i = 0, … n) p(xi)=yi (i=0,…n), 称为 拉格朗日插值多项式,记作 L n ( x ) L_n(x) Ln(x),即
L n ( x ) = ∑ j = 0 n y j l j ( x ) = ∑ j = 0 n y j ∏ i = 0 i ≠ j n x − x i x j − x i L_{n}(x)=\sum_{j=0}^{n} y_{j} l_{j}(x)=\sum_{j=0}^{n} y_{j} \prod_{{i=0 \\ i \neq j}}^{n} \frac{x-x_{i}}{x_{j}-x_{i}} Ln(x)=j=0∑nyjlj(x)=j=0∑nyji=0i=j∏nxj−xix−xi
三次样条插值
样条函数: 由一些按照某种光滑条件分段拼接起来的多项式组成的函数, 最常用的样条函数为三次样条函数,即由三次多项式组成,满足处处有二阶连续导数.
节点:
x
0
<
x
1
<
⋯
<
x
n
−
1
<
x
n
x_0 < x_1 < \cdots < x_n-1 < x_n
x0<x1<⋯<xn−1<xn,
函数值:
y
i
,
i
=
0
,
1
,
2
,
⋯
,
n
y_i, \ i = 0, 1, 2, \cdots , n
yi, i=0,1,2,⋯,n
确定三次样条函数 s ( x ) s(x) s(x) 满足 s ( x i ) = y i ( i = 0 , 1 , 2 , … , n ) s\left(x_{i}\right)=y_{i} \ (i=0,1,2, \ldots, n) s(xi)=yi (i=0,1,2,…,n)
由定义可设
s ( x ) = { s 1 ( x ) , x ∈ [ x 0 , x 1 ] s 2 ( x ) , x ∈ [ x 1 , x 2 ] ⋮ s n ( x ) , x ∈ [ x n − 1 , x n ] s(x)=\left\{\begin{array}{cl} s_{1}(x), & x \in\left[x_{0}, x_{1}\right] \\ s_{2}(x), & x \in\left[x_{1}, x_{2}\right] \\ \vdots & \\ s_{n}(x), & x \in\left[x_{n-1}, x_{n}\right] \end{array}\right. s(x)=⎩⎪⎪⎪⎨⎪⎪⎪⎧s1(x),s2(x),⋮sn(x),x∈[x0,x1]x∈[x1,x2]x∈[xn−1,xn]
其中
s
k
(
x
)
s_{k}(x)
sk(x) 为
[
x
k
−
1
,
x
k
]
\left[x_{k-1}, x_{k}\right]
[xk−1,xk] 上的三次多项式, 且满足
s
k
(
x
k
−
1
)
=
y
k
−
1
,
s
k
(
x
k
)
=
y
k
(
k
=
1
,
2
,
…
,
n
)
s_{k}\left(x_{k-1}\right)=y_{k-1}, s_{k}\left(x_{k}\right)=y_{k} \quad(k=1,2, \ldots, n)
sk(xk−1)=yk−1,sk(xk)=yk(k=1,2,…,n)
三次样条函数的连接条件:
s ( x ) ∈ C 2 [ a , b ] s(x) \in C^{2}[a, b] s(x)∈C2[a,b]
⟹ s ′ ( x k − ) = s ′ ( x k + ) , s ′ ′ ( x k − ) = s ′ ′ ( x k + ) \Longrightarrow s^{\prime}\left(x_{k}^{-}\right)=s^{\prime}\left(x_{k}^{+}\right), s^{\prime \prime}\left(x_{k}^{-}\right)=s^{\prime \prime}\left(x_{k}^{+}\right) ⟹s′(xk−)=s′(xk+),s′′(xk−)=s′′(xk+)
⟹ s k ′ ( x k − ) = s k + 1 ′ ( x k + ) , s k ′ ′ ( x k − ) = s k + 1 ′ ′ ( x k + ) \Longrightarrow s_{k}^{\prime}\left(x_{k}^{-}\right)=s_{k+1}^{\prime}\left(x_{k}^{+}\right), s_{k}^{\prime \prime}\left(x_{k}^{-}\right)=s_{k+1}^{\prime \prime}\left(x_{k}^{+}\right) ⟹sk′(xk−)=sk+1′(xk+),sk′′(xk−)=sk+1′′(xk+)
k = 1 , 2 , … , n − 1 k=1,2, \ldots, n-1 k=1,2,…,n−1
每个 s k ( x ) s_k(x) sk(x) 均为三次多项式,有 4 4 4 个待定系数,所以共有 4 n 4_n 4n 个待定系数,需 4 n 4_n 4n 个方程才能确定。前面已经得到 2 n + 2 ( n − 1 ) = 4 n – 2 2n +2(n -1) = 4n –2 2n+2(n−1)=4n–2 个方程,还缺 2 2 2 个方程,实际问题通常对样条函数在端点处的状态有要求,即所谓的边界条件
三次样条函数的边界条件:
- 第一类边界条件:给定函数在端点处的一阶导数,即
s ′ ( x 0 ) = f 0 ′ , s ′ ( x n ) = f n ′ s^{\prime}\left(x_{0}\right)=f_{0}^{\prime}, s^{\prime}\left(x_{n}\right)=f_{n}^{\prime} s′(x0)=f0′,s′(xn)=fn′
- 第二类边界条件:给定函数在端点处的二阶导数,即
s ′ ′ ( x 0 ) = f 0 ′ ′ , s ′ ′ ( x n ) = f n ′ ′ s^{\prime \prime}\left(x_{0}\right)=f_{0}^{\prime \prime}, s^{\prime \prime}\left(x_{n}\right)=f_{n}^{\prime \prime} s′′(x0)=f0′′,s′′(xn)=fn′′
当 s ′ ′ ( x 0 ) = s ′ ′ ( x n ) = 0 s''({x_0}) = s''({x_n}) = 0 s′′(x0)=s′′(xn)=0 时, 称为自然边界条件,此时的样条函数称为自然样条函数
- 第三类边界条件:设 f ( x ) f (x) f(x) 是周期函数,并设 x n – x 0 x_n – x_0 xn–x0 是一个周期,于是 s ( x ) s(x) s(x) 满足
s ′ ( x 0 ) = s ′ ( x n ) , s ′ ′ ( x 0 ) = s ′ ′ ( x n ) s^{\prime}\left(x_{0}\right)=s^{\prime}\left(x_{n}\right), s^{\prime \prime}\left(x_{0}\right)=s^{\prime \prime}\left(x_{n}\right) s′(x0)=s′(xn),s′′(x0)=s′′(xn)
拟合
- 插值: 过点 (适合精确数据)
- 拟合: 不过点, 整体近似;(适合有经验公式或有误差的数据)
直线最小二乘拟合
对于给定的数据点 ( x i , y i ) \left(x_{i}, y_{i}\right) (xi,yi), i = 1 , 2 , ⋯ , N i=1,2, \cdots, N i=1,2,⋯,N 求 y = a + b x y=a+b x y=a+bx, 使总误差平方和为最小, 即
Q = ∑ i = 1 N [ y i − ( a + b x i ) ] 2 Q=\sum_{i=1}^{N}\left[y_{i}-\left(a+b x_{i}\right)\right]^{2} Q=i=1∑N[yi−(a+bxi)]2
最小. Q Q Q 是关于未知数 a a a 和 b b b 的二元函数, 这一问题就是要确定 a a a 和 b b b 取何值时,二元函数 Q ( a , b ) Q(a,b) Q(a,b) 的值最小.
∂ Q ∂ a = 0 , ∂ Q ∂ b = 0 \frac{\partial \boldsymbol{Q}}{\partial a}=\mathbf{0}, \quad \frac{\partial \boldsymbol{Q}}{\partial b}=\mathbf{0} ∂a∂Q=0,∂b∂Q=0
∂ Q ∂ a = ∑ i = 1 N [ y i − ( a + b x i ) ] ⋅ ( − 1 ) = 0 ⟹ N a + b ∑ i = 1 N x i = ∑ i = 1 N y i \frac{\partial Q}{\partial a}=\sum_{i=1}^{N}\left[y_{i}-\left(a+b x_{i}\right)\right] \cdot(-1)=0 \Longrightarrow N a+b \sum_{i=1}^{N} x_{i}=\sum_{i=1}^{N} y_{i} ∂a∂Q=i=1∑N[yi−(a+bxi)]⋅(−1)=0⟹Na+bi=1∑Nxi=i=1∑Nyi
∂ Q ∂ b = ∑ i = 1 N [ y i − ( a + b x i ) ] ⋅ ( − x i ) = 0 ⟹ a ∑ i = 1 N x i + b ∑ i = 1 N x i 2 = ∑ i = 1 N x i y i \frac{\partial Q}{\partial b}=\sum_{i=1}^{N}\left[y_{i}-\left(a+b x_{i}\right)\right] \cdot\left(-x_{i}\right)=0 \Longrightarrow a \sum_{i=1}^{N} x_{i}+b \sum_{i=1}^{N} x_{i}^{2}=\sum_{i=1}^{N} x_{i} y_{i} ∂b∂Q=i=1∑N[yi−(a+bxi)]⋅(−xi)=0⟹ai=1∑Nxi+bi=1∑Nxi2=i=1∑Nxiyi
从而,未知数 a , b a, b a,b 满足如下方程组 (正则方程组)
{ N a + b ∑ i = 1 N x i = ∑ i = 1 N y i a ∑ i = 1 N x i + b ∑ i = 1 N x i 2 = ∑ i = 1 N x i y i \begin{cases} N a+b \sum_{i=1}^{N} x_{i}=\sum_{i=1}^{N} y_{i} \\ a \sum_{i=1}^{N} x_{i}+b \sum_{i=1}^{N} x_{i}^{2}=\sum_{i=1}^{N} x_{i} y_{i} \end{cases} {Na+b∑i=1Nxi=∑i=1Nyia∑i=1Nxi+b∑i=1Nxi2=∑i=1Nxiyi
解出 a , b a, b a,b, 就可得最小二乘拟合直线 y = a + b x y=a+bx y=a+bx
一般线性最小二乘
给定一个线性无关的函数系 ${ {\varphi _k}(x) \ | \ k = 1,2, \cdots ,m} $ ,如果拟合函数以其线性组合的形式
f ( x ) = ∑ k = 1 m a k φ k ( x ) f(x) = \sum\limits_{k = 1}^m {{a_k}{\varphi _k}(x)} f(x)=k=1∑makφk(x)
出现,则 f ( x ) = f ( x , a 1 , a 2 , ⋯ , a m ) f(x) = f(x,{a_1},{a_2}, \cdots ,{a_m}) f(x)=f(x,a1,a2,⋯,am) 就是关于参数 a 1 , a 2 , ⋯ , a m {a_1},{a_2}, \cdots ,{a_m} a1,a2,⋯,am 的线性函数,代入
J = ∑ i = 1 n ( f ( x i ) − y i ) 2 J = \sum\limits_{i = 1}^n {{{(f({x_i}) - {y_i})}^2}} J=i=1∑n(f(xi)−yi)2
则函数 J = J ( a 1 , a 2 , ⋯ , a k ) J = J({a_1},{a_2}, \cdots ,{a_k}) J=J(a1,a2,⋯,ak) 是关于参数 a 1 , a 2 , ⋯ , a m {a_1},{a_2}, \cdots ,{a_m} a1,a2,⋯,am 的多元函数。由
∂ J ∂ a k = 0 , k = 1 , 2 , ⋯ , m \frac{{\partial J}}{{\partial {a_k}}} = 0,\quad k = 1,2, \cdots ,m ∂ak∂J=0,k=1,2,⋯,m
亦即
∑ i = 1 n [ ( f ( x i ) − y i ) φ k ( x i ) ] = 0 , k = 1 , 2 , ⋯ , m \sum\limits_{i = 1}^n {[(f({x_i}) - {y_i}){\varphi _k}({x_i})]} = 0,\quad k = 1,2, \cdots ,m i=1∑n[(f(xi)−yi)φk(xi)]=0,k=1,2,⋯,m
可得
∑ j = 1 m [ ∑ i = 1 n φ j ( x i ) φ k ( x i ) ] a j = ∑ i = 1 n y i φ k ( x i ) , k = 1 , 2 , ⋯ , m \sum\limits_{j = 1}^m {\left[ {\sum\limits_{i = 1}^n {{\varphi _j}({x_i}){\varphi _k}({x_i})} } \right]{a_j}} = \sum\limits_{i = 1}^n {{y_i}{\varphi _k}({x_i})} ,\quad k = 1,2, \cdots ,m j=1∑m[i=1∑nφj(xi)φk(xi)]aj=i=1∑nyiφk(xi),k=1,2,⋯,m
记
R = ( φ 1 ( x 1 ) φ 2 ( x 1 ) ⋯ φ m ( x 1 ) φ 1 ( x 2 ) φ 2 ( x 2 ) ⋯ φ m ( x 2 ) ⋮ ⋮ ⋮ φ 1 ( x n ) φ 2 ( x n ) ⋯ φ m ( x n ) ) {\mathbf{R}} = \left( {\begin{array}{l} {{\varphi _1}({x_1})}&{{\varphi _2}({x_1})}& \cdots &{{\varphi _m}({x_1})} \\ {{\varphi _1}({x_2})}&{{\varphi _2}({x_2})}& \cdots &{{\varphi _m}({x_2})} \\ \vdots & \vdots &{}& \vdots \\ {{\varphi _1}({x_n})}&{{\varphi _2}({x_n})}& \cdots &{{\varphi _m}({x_n})} \end{array}} \right) R=⎝⎜⎜⎜⎛φ1(x1)φ1(x2)⋮φ1(xn)φ2(x1)φ2(x2)⋮φ2(xn)⋯⋯⋯φm(x1)φm(x2)⋮φm(xn)⎠⎟⎟⎟⎞
A = ( a 1 a 2 ⋮ a m ) {\mathbf{A}} = \left( {\begin{array}{c} a_1 \\ {{a_2}} \\ \vdots \\ {{a_m}} \end{array}} \right) A=⎝⎜⎜⎜⎛a1a2⋮am⎠⎟⎟⎟⎞
Y = ( y 1 y 2 ⋮ y n ) {\mathbf{Y}} = \left( {\begin{array}{c} {{y_1}} \\ {{y_2}} \\ \vdots \\ {{y_n}} \end{array}} \right) Y=⎝⎜⎜⎜⎛y1y2⋮yn⎠⎟⎟⎟⎞
则正规方程组可表示为
R
T
R
A
=
R
T
Y
{{\mathbf{R}}^T}{\mathbf{RA}} = {{\mathbf{R}}^T}{\mathbf{Y}}
RTRA=RTY
由代数知识可知,当矩阵
R
{\mathbf{R}}
R 是列满秩时
R
T
R
{{\mathbf{R}}^T}{\mathbf{R}}
RTR 是可逆的。于是正规方程组有唯一解,即
A
=
(
R
T
R
)
−
1
R
T
Y
{\mathbf{A}} = {({{\mathbf{R}}^T}{\mathbf{R}})^{ - 1}}{{\mathbf{R}}^T}{\mathbf{Y}}
A=(RTR)−1RTY
为所求的拟合函数的系数,就可得到最小二乘拟合函数 f ( x ) f(x) f(x)