2021-03-06-数值分析-Day09-分段插值-三次样条函数

5.6 - 5.7 分段插值

Runge现象：在一组等间插值点上使用具有高次多项式的多项式插值时出现的区间边缘处的振荡问题。

分段Lagrange插值

分段线性插值

通过相邻两个插值点作线性插值。已知节点 \(a=x_0<x_1<\dots<x_n=b\) ，记 \(h_k=x_{k+1}-x_k,\ h=\max\limits_{0\le k\le n-1}h_k\) ，记分段插值函数为 \(I_h(x)\) ，为 \(n-1\) 段折线。

余项估计有 \(|f(x)-I_h(x)|\le \displaystyle\frac{\max\limits_{a\le x\le b}|f''(x)|}{8}h^2\) ，说明分段线性插值函数具有一致收敛性。

分段二次插值

相邻三个插值点作二次插值。

课程中以线性插值中两邻点中间值，作为补充条件，即取 \(x_{i-0.5}=(x_i+x_{i-1})/2\) 作为每组中的第三个插值点，使得两点间距为 \(h_k/2\) 。按照这种方式得到余项估计 \(|f(x)-I_h(x)|\le \displaystyle\frac{\max\limits_{a\le x\le b}|f'''(x)|}{72\sqrt{3}}h^3\) 。

PS：个人认为应直接取连续三个插值点，即 \(\{x_0,x_1,x_2\}\ \{x_2,x_3,x_4\}\ \{x_4,x_5,x_6\}\ \dots\) 作为插值点进行二次插值。

分段Lagrange插值的问题：区间内出现不可导点

分段Hermite插值

设节点 \(a\le x_0<x_1<\dots<x_n\le b,\ h_i=x_i-x_{i-1}\ (i=1,2,\dots,n)\) ，给出插值条件： \(y_k=f(x_k),\ y'_k=f'(x_k),\ (k=0,1,\dots,n)\) 。则每区间 \([x_{i-1},x_i]\) 具有四个插值条件。构造三次多项式 \(H_3^{(i)}(x)\) ： \[ \begin{array}{c}H_3^{(i)}(x)=\varphi_{i-1}(x)y_{i-1}+\varphi_i(x)y_i+\psi_{i-1}(x)y'_{i-1}+\psi_i(x)y'_i\\ \left\{\begin{aligned} &\varphi_{i-1}(x_{i-1})=1 &&\varphi_{i-1}(x_i)=0 &&\varphi'_{i-1}(x_{i-1})=0 &&\varphi'_{i-1}(x_i)=0\\ &\varphi_i(x_{i-1})=0 &&\varphi_i(x_i)=1 &&\varphi'_i(x_{i-1})=0 &&\varphi'_i(x_i)=0\\ &\psi_{i-1}(x_{i-1})=0 &&\psi_{i-1}(x_i)=0 &&\psi'_{i-1}(x_{i-1})=1 &&\psi'_{i-1}(x_i)=0\\ &\psi_i(x_{i-1})=0 &&\psi_i(x_i)=0 &&\psi'_i(x_{i-1})=0 &&\psi'_i(x_i)=1\\ \end{aligned}\right.\\\Longrightarrow\quad\left\{\begin{aligned} &H_3^{(i)}(x_{i-1})=y_{i-1}&&H_3^{(i)}(x_i)=y_i\\ &H_3'^{(i)}(x_{i-1})=y'_{i-1}&&H_3'^{(i)}(x_i)=y'_i \end{aligned}\right.\end{array} \] 其中，称 \(\varphi_{i-1}(x),\ \varphi_i(x),\ \psi_{i-1}(x),\ \psi_i(x)\) 为三次Hermite插值基函数。
求得： \[ \begin{array}{c}\left\{\begin{aligned}\displaystyle &\varphi_{i-1}(x_{i-1})=\frac{1}{h_i^3}(2x-3x_{i-1}+x_i)(x-x_i)^2\\ &\varphi_i(x_{i-1})=\frac{1}{h_i^3}(-2x-x_{i-1}+3x_i)(x-x_{i-1})^2\\ &\psi_{i-1}(x_{i-1})=\frac{1}{h_i^2}(x-x_{i-1})(x-x_i)^2\\ &\psi_i(x_{i-1})=\frac{1}{h_i^2}(x-x_{i-1})^2(x-x_i) \end{aligned}\right.\\ \begin{aligned}\displaystyle H_3^{(i)}(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)y'_{i-1}+(x-x_{i-1})y'_i] }{h_i^2} \end{aligned}\end{array} \]

余项分析：若 \(f(x)\in C^4[a,b]\) （四阶连续可微），则当 \(x\in[x_{i-1},x_i]\) 时，有： \[ \begin{array}{c}\displaystyle f(x)-H_3^{(i)}(x)=\frac{f^{(4)}(\xi_i)}{4!}(x-x_{i-1})^2(x-x_i)^2,\quad\xi\in[x_{i-1},x_i]\\ \displaystyle |f(x)-H_3(x)|\le\frac{\max\limits_{a\le x\le b}|f^{(4)}(x)|\times h^4}{4!\times16}=\frac{\max\limits_{a\le x\le b}|f^{(4)}(x)|\times h^4}{384}\end{array} \] 可知 \(H_3(x)\) 是收敛的，且在 \([a,b]\) 内具有一阶连续导数。

5.8 - 5.10 三次样条插值

Hermite分段插值只能保证一阶连续可导，引入三次样条插值保证二阶连续可导。

三次样条函数：若函数 \(S(x)\in C^2[a,b]\) ，且在每个小区间 \([x_j,x_{j+1}]\) 上是三次多项式，其中 \(a=x_0<x_1<\dots<x_n=b\) 是给定节点，则称 \(S(x)\) 是节点 \(x_0,x_1,\dots,x_n\) 上的三次样条函数。

分段三次多项式： \(S(x)=a_jx^3+b_jx^2+c_jx+d_j,\ j=0,1,2,\dots,n-1\) ，其中 \(a_j\ b_j\ c_j\ d_j\) 为待定系数，因此函数共有 \(4n\) 个待定参数。函数满足关系： \(\left\{\begin{array}{l}S(x_j-0)=S(x_j+0)\\S'(x_j-0)=S'(x_j+0)\\S''(x_j-0)=S''(x_j+0)\end{array}\right.\) ，共 \(3n-3\) 个条件，加上插值条件 \(S(x_j)=f_j,\ j=0,1,\dots,n\) 共 \(4n-2\) 个条件。在此基础上在区间端点处添加一对边界条件即可满足唯一确定 \(4n\) 个待定参数的 \(4n\) 个条件。

边界条件：一般来说有两种边界条件，即两端点的一阶导数值或二阶导数值。

周期边界条件：当 \(f(x)\) 是以 \(x_n-x_0\) 为周期的周期函数时，要求 \(S(x)\) 也是周期函数，此时边界条件满足： \(\left\{\begin{array}{l}S(x_0+0)=S(x_n-0)\\S'(x_0+0)=S'(x_n-0)\\S''(x_0+0)=S''(x_n-0)\end{array}\right.\) ，且端点二阶导数值均为 \(0\) 。称这样的样条函数为周期样条函数。

代码：Hermite插值

//
// Created by xa on 2021-03-06.
//

#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>

using std::vector;

double x_im1; double x_i;
double y_im1; double y_i;
double y_1_im1; double y_1_i;

// x, y, y_1
vector<vector<double>> array;

double H(double x);

int main() {
    vector<vector<double>> array_ {{100, 10, 1/20}, {121, 11, 1/22}, {144, 12, 1/24}, {169, 13, 1/26}};
    array.assign(array_.begin(), array_.end());
    std::cout << H(125);
}

double H(double x) {
    sort(array.begin(), array.end(), [](const vector<double> &a, const vector<double> &b) {return a[0] < b[0]; });
    for (int i = 0; i < array.size(); ++i) {
        if ( x < array[i][0]) continue;
        else {
            x_im1 = array[i][0]; x_i = array[i+1][0];
            y_im1 = array[i][1]; y_i = array[i+1][1];
            y_1_im1 = array[i][2]; y_1_i = array[i+1][2];
            break;
        }
    }
    double h_i = x_i - x_im1;
    double h = ( ( (2 * x - 3 * x_im1 + x_i) * pow((x - x_i), 2) * y_im1
                 + (- 2 * x - x_im1 + 3 * x_i) * pow((x - x_im1), 2) * y_i )
             + h_i * (x - x_im1) * (x - x_i) * ((x - x_i) * y_1_im1 + (x - x_im1) * y_1_i ) )
             / pow(h_i, 3);
    return h;
}