续: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\)

代码:本节计算主要过程为解线性方程组,列出对应元计算式,调用对应函数即可。线性方程组解法在之前章节已给出,不再重复。