保形三次hermit插值

一、算法实现

一、插值函数建立

设函数 $y=F(x)$在区间$[a,b]$上有定义，且已知在离散点$a=x_0<x_1<...<x_n=b$上的值$y_0,y_1,...y_n，$$f(x)$在$[x_j,x_{j+1}]$分段区间内可表示为
$$f(x)=a(x-x_j{)}^3+b(x-x_j{)}^2+c(x-x_j)+d$$
设 $f^{\prime }(x)$是一阶导数，则
$$f^{\prime }(x)=3a(x-x_j{)}^2+2b(x-x_j)+c$$
将端点处 $f(x_j)=y_j$，$f(x_{j+1})=y_{j+1}$ 带入得
$$\{\begin{array}{rl}f(x_j)& =d\\ f^{\prime }(x_j)& =c\\ f(x_{j+1})& =a(x_{j+1}-x_j{)}^3+b(x_{j+1}-x_j{)}^2+c(x_{j+1}-x_j)\\ f^{\prime }(x_{j+1})& =3a(x_{j+1}-x_j{)}^2+2b(x_{j+1}-x_j)+c\end{array}$$
可得关于$a,b$ 的方程为
$$\{\begin{array}{rl}a(x_{j+1}-x_j{)}^3+b(x_{j+1}-x_j{)}^2& =f(x_{j+1})-f(x_j)-f^{\prime }(x_j)-f^{\prime }(x_j)(x_{j+1}-x_j)\\ 3a(x_{j+1}-x_j{)}^2+2b(x_{j+1}-x_j)& =f^{\prime }(x_{j+1})-f^{\prime }(x_j)\end{array}$$
记 $x_j$ 处的差商 ${\delta }_j=\frac{f(x_{j+1})-f(x_j)}{x_{j+1}-x_j}$，$x_j$ 处的一阶导 $d_j=f^{\prime }(x_j)$，$dzzdx=\frac{{\delta }_j-d_j}{x_{j+1}-x_j}，dzdxdx=\frac{d_{j+1}-{\delta }_j}{x_{j+1}-x_j}$，上式可表示为
$$\{\begin{array}{rl}a(x_{j+1}-x_j)+b& =dzzdx\\ 3a(x_{j+1}-x_j)+2b& =dzdxdx+dzzdx\end{array}$$
解方程组得
$$\{\begin{array}{rl}a& =\frac{dzdxdx-dzzdx}{x_{j+1}-x_j}\\ b& =2dzzdx-dzdxdx\end{array}$$

令 $h_j=x_{j{+}_1}-x_j$，最终得出
$$\{\begin{array}{rl}a& =\frac{d_{j+1}+d_j-2{\delta }_j}{{h_j}^2}\\ b& =\frac{-d_{j+1}-2d_j+3{\delta }_j}{h_j}\\ c& =f^{\prime }(x_j)=d_j\\ d& =f(x_j)=y_j\end{array}$$
其中 $h_j、{\delta }_j、y_j$ 均已知，求出 $x_j、x_{j+1}$ 处的导数 $d_j、d_{j+1}$ 方程得解

二、一阶导数求法

一、内点处的导数求法

内点处的一阶导数有以下规则：

1. 如果第 $k$ 个节点附近的差商 ${\delta }_{k-1}$ 和 ${\delta }_k$ 符号相反，或者其中一个为0，则该点处的一阶导数 $d_k=0$
2. 如果第 $k$ 个节点附近的差商 ${\delta }_{k-1}$ 和 ${\delta }_k$ 符号相同，则改点处的导数
   $$d_{k+1}=\frac{\delta mi{n}_k}{w{1}_k\frac{{\delta }_k}{\delta ma{x}_k}+w{2}_k\frac{{\delta }_{k+1}}{\delta ma{x}_k}}$$
   $其中h_k=(x_{k+1}-x_k)，h{s}_k=h_k+h_{k+1}，{\delta }_k=\frac{y_{k+1}-y_k}{x_{k+1}-x_k}，\delta mi{n}_k=min({\delta }_k,{\delta }_{k+1})，\delta ma{x}_k=max({\delta }_k,{\delta }_{k+1})，w{1}_k=\frac{h_k+h{s}_k}{3h{s}_k}，w{2}_k=\frac{h_{k+1}+h{s}_k}{3h{s}_k}$

二、端点处的导数求法

$$\{\begin{array}{rl}{d}_0& =\frac{(2h_0+h_1){\delta }_0-h_0{\delta }_1}{(h_0+h_1)}\\ d_n& =\frac{(2h_{n-1}+h_{n-2}){\delta }_{n-1}-h_{n-1}{\delta }_{n-2}}{(h_{n-2}+h_{n-1})}\end{array}$$

二、实验仿真

# -*- encoding: utf-8 -*-
'''
@File    :   pchip.py
@Time    :   2023/03/01 11:40:41
@Author  :   answer
'''# here put the import lib
import numpy as np
from matplotlib import pyplot as plt
from scipy import interpolatedef find_0point(delta):k = []for i in range(len(delta)-1):if delta[i] * delta[i+1] > 0:k.append(i)return k# 三次分段hermit函数
def pchip_spline(x, y, frequence):# x,y差分x_diff = []y_diff = []delta = []for i in range(len(x)-1):x_diff.append(x[i+1] - x[i])y_diff.append(y[i+1] - y[i])delta.append(y_diff[i]/x_diff[i])# 节点导数n = len(x)slope = [0 for i in range(n)]if n == 2:slope = [delta[0] for i in range(n)]else:k = find_0point(delta)for i in range(len(k)):index = k[i]dx_diff = x_diff[index] + x_diff[index + 1]w1 = (x_diff[index] + dx_diff) / (3 * dx_diff)w2 = (x_diff[index + 1] + dx_diff) / (3 * dx_diff)dmax = max(abs(delta[index]), abs(delta[index+1]))dmin = min(abs(delta[index]), abs(delta[index+1]))slope[index + 1] = dmin / \(w1*delta[index]/dmax + w2*delta[index+1]/dmax)slope[0] = 0# 库函数默认端点导数不为0 interpolate.pchip_interpolate(x, y, x_pchip)# slope[0] = ((2 * x_diff[0] + x_diff[1]) * delta[0] -#             x_diff[0] * delta[1]) / (x_diff[0] + x_diff[1])# if slope[0] * delta[0] < 0:#     slope[0] = 0# elif (delta[0] * delta[1] < 0) & (abs(slope[0]) > 3 * abs(delta[0])):#     slope[0] = 3 * delta[0]# print(slope)# slope[n - 1] = ((2 * x_diff[n - 2] + x_diff[n - 3]) * delta[n - 2] -#                 x_diff[n - 2] * delta[n - 3]) / (x_diff[n - 3] + x_diff[n - 2])# if delta[n - 2] * slope[n - 1] < 0:#     slope[n - 1] = 0# elif (delta[n - 2] * delta[n - 3] < 0) & (abs(slope[n - 1]) > 3 * abs(delta[n - 2])):#     slope[n - 1] = 3 * delta[n - 2]# print(slope)# hermit splinex_hermit = []y_hermit = []for i in range(n - 1):# 计算多项式系数a = (slope[i + 1] + slope[i] - 2 * delta[i]) / (x_diff[i]**2)b = (3 * delta[i] - 2 * slope[i] - slope[i + 1]) / x_diff[i]c = slope[i]d = y[i]# 计算插值点for j in range(frequence):x_inter = x[i] + j * (x[i+1] - x[i]) / frequencex_hermit.append(x_inter)y_hermit.append(a * (x_inter - x[i])**3 + b * (x_inter - x[i])**2 + c * (x_inter - x[i]) + d)x_hermit.append(x[n-1])y_hermit.append(y[n-1])return x_hermit, y_hermitif __name__ == '__main__':frequence = 10x = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]y = [0, 1.5, 0, 0, 0.5, 0.4, 1.2, 1.2,0.1, 0, 0.3, 0.6]x_pchip, y_pchip = pchip_spline(x, y, frequence)y_ = interpolate.pchip_interpolate(x, y, x_pchip)y_1 = interpolate.splrep(x, y)y_1 = interpolate.splev(x_pchip, y_1)y_2 = interpolate.Akima1DInterpolator(x, y)y_2 = y_2(x_pchip)y_3 = interpolate.interp1d(x, y, 'cubic')y_3 = y_3(x_pchip)plt.subplot(2, 2, 1)plt.plot(x, y, "o", label='sample point')plt.plot(x_pchip, y_, linewidth=3.0, label="scipy_pchip")plt.plot(x_pchip, y_1, color='lime', label="spline")plt.legend()plt.subplot(2, 2, 2)plt.plot(x, y, "o", label='sample point')plt.plot(x_pchip, y_, linewidth=3.0, label="scipy_pchip")plt.plot(x_pchip, y_3, color='blueviolet', label="cubic")plt.legend()plt.subplot(2, 2, 3)plt.plot(x, y, "o", label='sample point')plt.plot(x_pchip, y_, linewidth=3.0, label="scipy_pchip")plt.plot(x_pchip, y_2, color='dodgerblue', label="akima")plt.legend()plt.subplot(2, 2, 4)plt.plot(x, y, "o", label='sample point')plt.plot(x_pchip, y_, linewidth=3.0, label="pchip_scipy")plt.plot(x_pchip, y_pchip, color='gray', label="pchip d_point=0")plt.legend()plt.subplots_adjust(left=0.04, bottom=0.05, right=0.98,top=0.98, wspace=0.08, hspace=0.1)plt.show()

1. 以 $(0,4)，(1,3)，(2,4)，(3,6)，(5,7)，(6,5)，(8,8)，(11,1)$ 作为节点，每点之间插十次
   

2. 对阶跃信号进行插值