续:5.8 - 5.10 三次样条插值
三转角方法
考虑第一种一般边界条件: \(S'(x_0)=f_0',\ S'(x_n)=f_n'\) ,即已知两端点一阶导数值。
令 \(m_i=S'(x_i),\
i=0,1,\dots,n\) ,利用三次Hermite插值,得到 \(S(x)=\displaystyle
\sum_{j=0}^n[y_j\alpha_j(x)+m_j\beta_j(x)]\) ,其中 \(\alpha_j(x),\ \beta_j(x)\)
为分段三次Hermite插值的基函数。再由边界条件得 \(S'(x_0)=f_0',\
S'(x_n)=f_n'\) 即可解出 \(m_i\) 在各插值点的取值。记
\(\displaystyle
\lambda_i=\frac{h_{i+1}}{h_i+h_{i+1}},\
\mu_i=1-\lambda_i=\frac{h_i}{h_i+h_{i+1}},\
g_i=3(\lambda_if[x_{i-1},x_i]+\mu_if[x_i,x_{i+1}])\) ,
最终解得: \[
\begin{array}{c}\lambda_im_{i-1}+2m_i+\mu_im_{i+1}=g_i\\
\begin{pmatrix}
2 &\mu_1 & & & & \\
\lambda_2 &2 &\mu_2 & & & \\
&\ddots &\ddots &\ddots & & \\
& &\ddots &\ddots &\ddots & \\
& & &\lambda_{n-2} &2 &\mu_{n-2} \\
& & & &\lambda_{n-1} &2
\end{pmatrix}
\begin{pmatrix}
m_1 \\
m_2 \\
\vdots \\
\vdots \\
m_{n-2} \\
m_{n-1}
\end{pmatrix}=
\begin{pmatrix}
g_1-\lambda_iy_0' \\
g_2 \\
\vdots \\
\vdots \\
g_{n-2} \\
g_{n-1}-\mu_{n-1}y_n'
\end{pmatrix}
\end{array}
\] 利用大型稀疏矩阵线性方程数值解法,解出 \(m_i\) ,即解得 \(x\in[x_{i-1},x_i]\) 时,有:
\(\begin{aligned}\displaystyle
S(x)&=\frac{
(2x-3x_{i-1}+x_i)(x-x_i)^2y_{i-1}+(-2x-x_{i-1}+3x_i)(x-x_{i-1})^2y_i
}{h_i^3}\\ &+\frac{
(x-x_{i-1})(x-x_i)[(x-x_i)m_{i-1}+(x-x_{i-1})m_i] }{h_i^2}
\end{aligned}\)
其他边界条件也可用类似方法解得,计算方法较复杂。
三弯矩方法
二阶导数 \(S''(x)=M_j\ (j=0,1,\dots,n)\) 在力学上解释为细梁在 \(x_j\) 截面处的弯矩。
令 \(M_i=S''(x_i),\
i=0,1,\dots,n\) ,则对 \(x\in[x_{i-1},x_i]\) ,有:
\(\displaystyle
S''(x)=\frac{x-x_i}{x_{i-1}-x_i}M_{i-1}+\frac{x-x_{i-1}}{x_i-x_{i-1}}M_i\)
,对此连续积分两次,记 \(S(x_i)=y_i\)
,得:
\(\begin{aligned}\displaystyle
S(x)&=\frac{1}{6h_i}\left[(x_i-x)^3M_{i-1}+(x-x_{i-1})^3M_i\right]\\&+\left(\frac{y_{i-1}}{h_i}-\frac{h_iM_{i-1}}{6}\right)(x_i-x)+\left(\frac{y_i}{h_i}-\frac{h_iM_i}{6}\right)(x-x_{i-1})\\&=\frac{(x_i-x)^3M_{i-1}+(x-x_{i-1})^3M_i+(6y_{i-1}-h_i^2M_{i-1})(x_i-x)+(6y_i-h_i^2M_i)(x-x_{i-1})}{6h_i}\end{aligned}\)
利用 \(S'(x_i-0)=S'(x_i+0)\)
求出 \(M_i\) ,记
\(\displaystyle
\lambda_i=\frac{h_{i+1}}{h_i+h_{i+1}}=,\ \mu_i=\frac{h_i}{h_i+h_{i+1}},\
d_i=6f[x_{i-1},x_i,x_{i+1}]\) ,
则有: \(\mu_iM_{i-1}+2M_i+\lambda_iM_{i+1}=d_i\)
结合边界条件 \(M_0=S''(x_0)=y''_0,\
M_n=S''(x_n)=y''_n\) ,可得: \[
\begin{array}{c}\mu_iM_{i-1}+2M_i+\lambda_iM_{i+1}=d_i\\
\begin{pmatrix}
2 &\mu_1 & & & & \\
\mu_2 &2 &\lambda_2 & & & \\
&\ddots &\ddots &\ddots & & \\
& &\ddots &\ddots &\ddots & \\
& & &\mu_{n-2} &2 &\lambda_{n-2} \\
& & & &\mu_{n-1} &2
\end{pmatrix}
\begin{pmatrix}
M_1 \\
M_2 \\
\vdots \\
\vdots \\
M_{n-2} \\
M_{n-1}
\end{pmatrix}=
\begin{pmatrix}
d_1-\mu_iy_0' \\
d_2 \\
\vdots \\
\vdots \\
d_{n-2} \\
d_{n-1}-\lambda_{n-1}y_n''
\end{pmatrix}
\end{array}
\] 或结合边界条件 \(S'(x_0)=y_0',\
S'(x_n)=y_n'\) ,得: \[
\begin{array}{c}\mu_iM_{i-1}+2M_i+\lambda_iM_{i+1}=d_i\\2M_0+M_1=d+0\\M_{n-1}+2M_n=d_n\\
\begin{pmatrix}
2 &_1 & & & & \\
\mu_1 &2 &\lambda_1 & & & \\
&\ddots &\ddots &\ddots & & \\
& &\ddots &\ddots &\ddots & \\
& & &\mu_{n-1} &2 &\lambda_{n-1} \\
& & & &1 &2
\end{pmatrix}
\begin{pmatrix}
M_1 \\
M_2 \\
\vdots \\
\vdots \\
M_{n-2} \\
M_{n-1}
\end{pmatrix}=
\begin{pmatrix}
d_0 \\
d_1 \\
\vdots \\
\vdots \\
d_{n-1} \\
d_n
\end{pmatrix}
\end{array}
\] 或结合周期边界条件 \(S'(x_0)=S'(x_n),\
S''(x_0)=S''(x_n)\) ,同样得到稀疏矩阵线性方程。
采用稀疏矩阵线性方程解法即可。
综上,当边界条件为第一类(已知一阶导数值)时,适合采用三转角方法;当边界条件为第二类(已知二阶导数值)时适合采用三弯矩方法;对周期边界条件,两种方法计算量一致。
最小二乘法
数据拟合问题
经由观察或测试得到 \(y(x)\) 的一组离散数据,在给定的函数类 \(\Phi\) 上根据这组离散数据做出逼近曲线,要求逼近曲线在 \(x_i\) 处与离散数据尽可能接近。
对函数 \(\varphi(x)\in\Phi\) ,要求以 \(\varphi(x)\) 在离散点的误差 \(\delta_i=\phi(x_i)-y_i,\ i=0,1,\dots,m\) 为分量的误差向量 \(\delta=(\delta_0,\delta_1,\dots,\delta_m)^T\) ,使某一向量范数 \(\|\delta\|\) 达到最小。对不同的范数,可以构造出不同意义下的拟合函数。
函数类 \(\Phi\) 通常取为 \(\Phi=Span\{\varphi_0(x),\varphi_1(x),\dots,\varphi_n(x)\}\) ,其中函数系 \(\varphi_0(x),\varphi_1(x),\dots,\varphi_n(x)\) 在包含节点 \(\{x_i\}\) 的区间 \([a,b]\) 上线性无关, \(\Phi\) 中任一函数 \(\varphi(x)\) 可表示为: \(\varphi(x)=a_0\varphi_0(x)+a_1\varphi_1(x)+\dots+a_n\varphi_n(x)\) 。
常用函数系包括幂函数系 \(\{x^j\}\) ,三角函数系 \(\{\sin jx\}\ \{\cos jx\}\) ,指数函数系 \(\{e^{\lambda_jx}\}\) ,正交函数系等。
在求误差向量 \(\delta\) 的范数时,常用 \(2\) -范数,对应的曲线拟合方法就称为最小二乘法。
推导:在函数类 \(\Phi\) 中找到函数 \(y=\varphi^*(x)\) ,使误差向量 \(\delta\) 的 \(2\)
-范数达到最小值。(通常取误差向量的加权形式 \(\rho(x)\varphi(x)\) 求 \(2\) -范数的平方 \(\|\delta\|_2^2\) )。该范数为关于 \((a_0,a_1,\dots,a_n)\) 的函数,记其为 \(G(a_0,a_1,\dots,a_n)\)
,则问题转换为求该函数最小值问题。
记 \(\left\{\begin{array}{l}\varphi_j=(\varphi_j(x_0),\varphi_j(x_1),\dots,\varphi_j(x_m))^T,\
j=0,1,2,\dots,n\\f=(y_0,y_1,\dots,y_m)^T\\(\varphi_j,\varphi_k)=\displaystyle\sum_{i=0}^m\rho(x_i)\varphi_j(x_i)\varphi_k(x_i)\\(f,\varphi_k)=\displaystyle\sum_{i=0}^m\rho(x_i)y_i\varphi_k(x_i)\end{array}\\\right.\)
得到关于系数向量 \((a_0,a_1,\dots,a_n)^T\) 的线性方程组:
\[
\begin{array}{c}\displaystyle\sum_{j=0}^n(\varphi_j,\varphi_k)a_j=(f,\varphi_k)\quad
k=0,1,\dots,n\\
\begin{pmatrix}
(\varphi_0,\varphi_0) &(\varphi_0,\varphi_1) &\cdots &(\varphi_j,\varphi_k)
\\
(\varphi_1,\varphi_0) &(\varphi_1,\varphi_1) &\cdots &(\varphi_1,\varphi_n)
\\
\vdots &\vdots &\ddots &\vdots \\
(\varphi_n,\varphi_0) &(\varphi_n,\varphi_1) &\cdots &(\varphi_n,\varphi_n)
\end{pmatrix}
\begin{pmatrix}
a_0 \\
a_1 \\
\vdots \\
a_n
\end{pmatrix}
=\begin{pmatrix}
(f,\varphi_0) \\
(f,\varphi_1) \\
\vdots \\
(f,\varphi_n)
\end{pmatrix}\end{array}
\]
称之为正则方程组或法方程组。
正则方程组的系数矩阵为对称矩阵。若向量组 \(\varphi_0,\varphi_1,\dots,\varphi_n\)
线性无关,则其为对称正定矩阵,可以通过迭代法求得唯一解 \(a_j^*\) ,从而得到拟合函数 \(\varphi^*(x)\) 。
取函数类
取函数类为幂函数系,则正则方程组为:
\[
\begin{pmatrix}
\sum\rho_i &\sum\rho_i x_i &\cdots &\sum\rho_i x_i^n \\
\sum\rho_i x_i &\sum\rho_i x_i^2 &\cdots &\sum\rho_i
x_i^{n+1} \\
\vdots &\vdots &\ddots &\vdots \\
\sum\rho_i x_i^n &\sum\rho_i
x_i^{n+1} &\cdots &\sum\rho_i x_i^{2n}
\end{pmatrix}
\begin{pmatrix}
a_0 \\
a_1 \\
\vdots \\
a_n
\end{pmatrix}
=\begin{pmatrix}
\sum\rho_iy_i \\
\sum\rho_ix_iy_i \\
\vdots \\
\sum\rho_ix_i^ny_i
\end{pmatrix}
\] 其中 \(\rho_i=\rho(x_i),\
\sum=\sum\limits_{i=0}^m\) 。此时拟合曲线 \(\varphi^*(x)=p_n^*(x)=a_0^*+a_1^*x+\dots+a_n^*x^n\)
,称之为多项式拟合曲线。
取函数类为正交函数系,则 \((\varphi_i,\varphi_j)=0,\ (i\neq j)\) ,正则方程组的系数矩阵为对角矩阵, \(a_k^*=(f,\varphi_k)/(\varphi_k,\varphi_k),\ k=0,1,\dots,n\) 。