插值是什么
在工程中，我们经常要用一条曲线将一些点依次连接起来，称为插值。
插值的可行性证明
插值法定理:对n+1个不同的节点有唯一多项式\(\phi (x)=a_0+a_1x+\cdots+a_nx^n\),使得\(\phi_n(x_j)=y_j(j=0,1,2,\cdots,n)\)
证明：将\(x_0\)到\(x_n\)带入多项式能得到一个线性方程组，AX=Y.其中A中元素为\(x_i\)，X中为\(a_i\)，Y中为\(y_i\).
A的形式是一个范德蒙行列式，一定可逆，所以线性方程有解，多项式一定存在。
插值函数
已知n+1个点\((x_i,y_i) (i=0,1,\cdots,n)\)下面讲各种插值函数。
分段线性插值函数
简单说就是将相邻的点两两用直线连接起来，如此形成一条折线，折线的方程为\(I_n(x)=\sum_{i=0}^ny_i l_i(x)\)，满足\(I_n(x_i)=y_i\)其中
\(l_i(x)=\left\{ \begin{matrix} \frac{x-x_{i-1}}{x_i-x_{i-1}},x\in [x_{i-1},x_i], i\neq 0 \\ \frac{x-x_{i+1}}{x_i-x_{i+1}},x\in [x_{i},x_{i+1}], i\neq n \\ 0, 其他 \end{matrix} \right.\)
相当于是在求每个点的y值时将前面的点累加起来。
拉格朗日插值
\(L_n(x)=\sum_{i=0}^ny_i l_i(x)=\sum_{i=0}^ny_i(\prod_{j=0,j\neq i}^n \frac{x-x_j}{x_i-x_j})\)
对于每一个\(x_i,都有l_i(x_i)\neq 0, l_j(x_i)=0 (j\neq i)\).
这样求\(I_n(x_i)=y_i\).
function y=lagrange(x0,y0,x)
%x0,y0为节点
%x是插值点
n = length(x0);m = length(x);
for i = 1:m
z = x(i);
s = 0.0;
for k = 1:n
p = 1.0;
for j = 1:n
if j ~= k
p = p * (z - x0(j)) / (x0(k) - x0(j));
end
end
s = p*y0(k)+s;
end
y(i)=s;
end
牛顿插值
拉格朗日插值法中存在的“增加一个节点时整个计算工作重新开始”的问题。
牛顿插值法是拉格朗日插值法的改进方法，运用迭代的思路求函数，且可以节省乘除法运算次数, 同时在Newton插值多项式中用到差分与差商等概念，又与数值计算的其他方面有密切的关系。
function [A,y]= newtonzi(X,Y,x)
% Newton插值函数
% X为已知数据点的x坐标
% Y为已知数据点的y坐标
% x为插值点的x坐标
% 函数返回A差商表
% y为各插值点函数值
n=length(X); m=length(x);
for t=1:m
z=x(t); A=zeros(n,n);A(:,1)=Y';
s=0.0; y=0.0; c1=1.0;
for j=2:n
for i=j:n
A(i,j)=(A(i,j-1)- A(i-1,j-1))/(X(i)-X(i-j+1));
end
end
C=A(n,n);
for k=1:n
p=1.0;
for j=1:k-1
p=p*(z-X(j));
end
s=s+A(k,k)*p;
end
ss(t)=s;
end
y=ss;
A=[X',A];
end
Hermite插值法
牛顿和拉格朗日插值法都有一些缺陷，它们只考虑了给出的点的拟合特性，没有考虑原函数的其他性质如导数与原函数是否适配。
有时候我们需要插值函数与原函数不仅需要一阶导数相同，还需要更高阶。
设原函数的导数为m(i),插值函数为H_i. 则Hermite函数满足
\(\left\{
\begin{matrix}
H(x_i)=y_i\\
H'(x_i)=m_i
\end{matrix}
\right.\)
具体算法见：https://blog.csdn.net/SanyHo/article/details/106849323
下面给出代码.
function f = Hermite( x,y,y_1,x0 )
%Hermite插值函数
% x为已知数据点的x坐标
% y为已知数据点的y坐标
% y_1为数据点y值导数
% x0为插值点的x坐标
syms t;
f = 0.0;
if(length(x) == length(y))
if(length(x) == length(y_1))
n = length(x);
else
disp('y和y的导数维数不相等');
renturn;
end
else
disp('x和y的维数不相等');
return;
end
%以上为输入判断和确定“n”的值
for i = 1:n
h = 1.0;
a = 0.0;
for j = 1:n
if(j ~= i)
h = h*(t-x(j))^2/((x(i)-x(j))^2);%求得值为(li(x))^2
a = a + 1/(x(i)-x(j)); %求得ai（x）表达式之中的累加部分
end
end
f = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));
if(i == n)
if(nargin == 4)
f = subs(f, 't' , x0); %输出结果
else
f = vpa(f,6); %输出精度为有效数字为6位的函数表达式
end
end
end
样条插值
给定区间[a,b]的一个划分
\(\Delta :a=x_0<x_1<\cdots <x_{n-1}<x_n=b\)
如果S(x)满足：
- 在每一个小区间上[\(x_i,x_{i+1}\)]，S(x)为m次多项式
- 它有m-1阶连续导数。
则S(x)为m次样条函数。
计算比较繁琐，就不详细说了，感兴趣看https://blog.csdn.net/qq_41035283/article/details/121459043
https://blog.csdn.net/qq_41035283/article/details/121459043
% 定义三次样条函数的系数矩阵、插值宽度、差商表、g值和M值
function [D,h,A,g,M]=threesimple(X,Y)
% X: 已知的插值节点数组
% Y: 插值节点对应的函数值数组
n = length(X); % 计算插值节点的数量
A = zeros(n,n); % 初始化差商表为n*n的零矩阵
A(:,1) = Y'; % 第一列填充插值点的函数值
D = zeros(n-2,n-2); % 初始化系数矩阵D为(n-2)*(n-2)的零矩阵
g = zeros(n-2,1); % 初始化g值为(n-2)*1的零向量
% 构建差商表A
for j = 2:n
for i = j:n
A(i,j) = (A(i,j-1) - A(i-1,j-1)) / (X(i) - X(i-j+1));
end
end
% 计算插值宽度h
for i = 1:n-1
h(i) = X(i+1) - X(i);
end
% 构建系数矩阵D和初始化g值
for i = 1:n-2
D(i,i) = 2;
g(i) = (6 / (h(i+1) + h(i))) * (A(i+2,2) - A(i+1,2));
end
% 填充D矩阵的非对角线元素
for i = 2:n-2
u(i) = h(i) / (h(i) + h(i+1));
n(i-1) = h(i) / (h(i-1) + h(i));
D(i-1,i) = n(i-1);
D(i,i-1) = u(i);
end
% 解线性方程组求解M值
M = D \\ g;
M = [0; M; 0]; % 边界条件，M的首尾元素设为0
end
% 定义三次样条插值函数
function s = threesimple1(X,Y,x)
% X: 已知的插值节点数组
% Y: 插值节点对应的函数值数组
% x: 需要计算插值值的点
[D,h,A,g,M] = threesimple(X,Y); % 调用threesimple函数计算相关参数
n = length(X); % 插值节点的数量
m = length(x); % 需要计算插值值的点的数量
% 计算每个插值点的函数值
for t = 1:m
for i = 1:n-1
if (x(t) <= X(i+1)) && (x(t) >= X(i))
% 计算三次样条插值公式的各个部分
p1 = M(i,1) * (X(i+1) - x(t)) ** 3 / (6 * h(i));
p2 = M(i+1,1) * (x(t) - X(i)) ** 3 / (6 * h(i));
p3 = (A(i,1) - M(i,1) / 6 * (h(i) ** 2)) * (X(i+1) - x(t)) / h(i);
p4 = (A(i+1,1) - M(i+1,1) / 6 * (h(i) ** 2)) * (x(t) - X(i)) / h(i);
s(t) = p1 + p2 + p3 + p4; % 计算插值点的函数值
break; % 找到对应的区间后跳出循环
else
s(t) = 0; % 如果x不在X的范围内，插值结果为0
end
end
end
end
matlab 工具包（香）
matlab有非常多的插值函数，详见官方文档：插值 - MATLAB & Simulink - MathWorks 中国
下面是一些例子。
三次Hermite插值法
p = pchip(x,y,new_x)
三次样条插值
p = spline(x,y,new_x)
n维插值
2维为例，函数为二元函数
p = interp2(x,y,z,new_x,new_y,method)
%method具体有 'linear','cublic','spline','nearest'
%x,y需要写成网格形式, [x,y] = meshgrid(x,y)