插值

插值：求过已知有限个数据点的近似函数

一维插值

插值多项式
已知函数$f(x)$在区间$[a,b]$上n+1个不同点$x_0,x_1,...,x_n$处的函数值$y_i=f(x_i)(i=0,1,...,n)$,
求一个至多n次多项式${\varphi }_n(x_i)=a_0+a_{1}x+...+a_n{x}^n$,使其满足${\varphi }_n(x_i)=f(x_i)(i=0,1,...,n)$，${\varphi }_n(x)$称为插值多项式
插值多项式与被插函数之间的差$R_n(x)=f(x)-{\varphi }_n(x)$，称为截断误差（插值余项）
f(x)充分光滑时有,$R_n(x)=f(x)-L_n(x)=\frac{f^{n+1}(\xi )}{(n+1)!}{w}_{n+1}(x),\xi \in (a,b)$,其中$w_{n+1}(x)={\prod }_{j=0}^n(x-x_j)$
解下列方程即可求出插值多项式

拉格朗日插值

1. 原理
   拉格朗日插值实际并未求出插值多项式系数，而是构造了一组基函数
   
2. matlab代码

%设n个节点数据以数组x0, y0输入
%m个插值点以数组x输入，输出数组y为m个插值
function y=lagrange(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

3. 示例
   

牛顿插值

原理
特点：每增加一个节点，插值多项式只增加一项

1. 差商
   (2) 一阶差商：$f[x_i,x_j]=\frac{f(x_i)-f(x_j)}{x_i-x_j}$
   (3) 二阶差商：$f[x_i,x_j,x_k]=\frac{f[x_i,x_j]-f[x_j,x_k]}{x_i-x_k}$
   (4) k阶差商：$f[x_0,x_1,...,x_k]=\frac{f[x_0,x_1,...,x_k-1]-f[x_1,x_2,...,x_k]}{x_0-x_k}$
2. Newton插值公式
   （1）线性公式 ：${\varphi }_1(x)=f(x_0)+(x-x_0)f[x_0,x_1]$
   （2）$N_n(x)=f(x_0)+(x-x_0)[x_0,x_1]+...+(x-x_0)(x-x_1)...(x-x_{n-1})f[x_0,x_1,...,x_n]$
   插值余项：$R_n(x)=(x-x_0)(x-x_1)...(x-x_n)f[x,x_0,x_1,...,x_n]=w_{n+1}(x)f[x,x_0,x_1,...,x_n]$
3. 差分(前向差分)
   结点等距(h)
   (1) 一阶差分：$\Delta f_k=f_{k+1}-f_k$
   (2) 二阶差分：${\Delta }^2{f}_k=\Delta f_{k+1}-\Delta f_k$
   (3) m阶差分：${\Delta }^m{f}_k={\Delta }^{m-1}{f}_{k+1}-{\Delta }^{m-1}{f}_k$
4. 等距结点插值公式
   牛顿前向插值公式：$N_n(x_0+th)=f_0+t\Delta f_0+...+\frac{t(t-1)...(t-n+1)}{n}{\Delta }^n{f}_0$, $(x=x_0+th)$

matlab代码

%newton.m
%求牛顿插值多项式、差商、插值及其误差估计的MATLAB主程序
%输入的量:X是n+1个节点(x_i,y_i)(i = 1,2, ... , n+1)横坐标向量，Y是纵坐标向量，
%x是以向量形式输入的m个插值点，M在[a,b]上满足｜f~(n+1)(x)｜≤M
%注：f~(n+1)(x)表示f(x)的n+1阶导数
%输出的量：向量y是向量x处的插值，误差限R，n次牛顿插值多项式L及其系数向量C，
%差商的矩阵A
function[y,R,A,C,L] = newton(X,Y,x,M)
n = length(X);
m = length(x);
for t = 1 : m
	z = x(t);
    	A = zeros(n,n);
	A(:,1) = Y';
    	s = 0.0; p = 1.0; q1 = 1.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
            q1 = abs(q1*(z-X(j-1)));
            c1 = c1 * j;
        end
        C = A(n, n); q1 = abs(q1*(z-X(n)));
        for k = (n-1): -1:1
            C = conv(C, poly(X(k)));
            d = length(C);
            C(d) = C(d) + A(k,k);%在最后一维，也就是常数项加上新的差商
        end
        y(t) = polyval(C,z);
        R(t) = M * q1 / c1;
end
L = poly2sym(C);

分段插值

1. 原理
   将每两个相邻的节点用直线连起来，形成的折线就是分段线性插值函数.

2. matlab实现
$y=interp1(x0,y0,x{,}^{\prime }linea{r}^{\prime })$

Hermite插值

1. 原理
   不仅要求节点处函数同值，而且要求与函数有相同的一阶、二阶、多阶导
   已知：$y=f(x)$在n+1个互异结点$x_0,x_1,...,x_n$上的函数值$y_i=f(x_i),(i=0,1,...,n)$和导数值$y_i^{\prime }=f^{\prime }(x_i),(i=0,1,...,n)$
   求：至多2n+1次的多项式$H(x)$，使得$H(x_i)=y_i,H^{\prime }(x_i)=y_i^{\prime },(i=0,1,...,n)$
   Hermite插值多项式：$H(x)={\sum }_{i=0}^n{h}_i[(x_i-x)(2a_i{y}_i-y_i^{\prime })+y_i]$
   其中$h_i={\prod }_{j=0,j\ne i}^n(\frac{x-x_j}{x_i-x_j}{)}^2,a_i={\sum }_{j=0,j\ne i}^n\frac{1}{x_i-x_j}$
2. matlab代码

n个节点的数据以数组x0（已知点的横坐标），y0（函数值），y1（导数值）输入
% m个插值点以数组x输入，输出数组y为m个插值
function y=hermite(x0,y0,y1,x)
n=length(x0);m=length(x); 
	for k=1:m 
 	yy=0.0; 
 		for i=1:n 
 			h=1.0; 
 			a=0.0; 
			for j=1:n 
 				if j~=i 
 					h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; 
 					a=1/(x0(i)-x0(j))+a; 
 				end 
 			end 
 			yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); 
 		end 
 	y(k)=yy; 
end

3. 示例
   

样条插值

对插值函数的光滑性有较高的要求
原理
给定区间[a,b]的一个划分$\Delta :a=x_0<x_1<...<x_{n-1}<x_n=b$
如果$s(x)$满足：
(i) 在每个小区间$[x_i,x_{i-1}](i=0,1,...,n-1)$上*s(x)*是$k$次多项式;
(ii) $s(x)$在$[a,b]$上具有$k-1$阶连续导数
则称$s(x)$为关于分划$\Delta$
的$k$次样条函数

1. 折线是一次样条函数
2. 二次样条函数
   
3. 三次样条函数
   
   matlab中的样条插值
   1. 三次样条曲线有4种边界条件:
   (i)自然边界条件，二阶导数在边界处为0，可视为简支梁，是最常用的边界条件。
   (ii)第一边界条件，二阶导在边界处已知的边界条件。自然边界条件可视为特例。
   (iii)第二边界条件，一阶导在边界处已知的边界条件。
   (iv)循环边界条件，一阶导与二阶导在边界处相等的边界条件，适用于封闭或循环的图形。
   (1) $y=interp1(x0,y0,x{,}^{\prime }splin{e}^{\prime })$
   (2) $y=spline(...)$
   spline函数只能实现自然边界条件和第二边界条件，可以实现一维或者高维的曲线插值
   (3) $pp=csape(x0,y0,conds)，y=ppval(pp,x)$
   csape函数可以实现三次样条曲线的各种条件，包括第一边界、第二边界、循环边界、混合边界等。可以实现一维至多维的各种情况。
   

matlab中的一维插值函数interp1

B 样条函数插值方法

二维插值

示 例：

clc,clear
[x,y]=meshgrid(-4:0.8:4); %原始数据
z=peaks(x,y);
[xi,yi]=meshgrid(-4:0.2:4) %插值数据
zi_nearest=interp2(x,y,z,xi,yi,‘nearest’);
zi_linear=interp2(x,y,z,xi,yi);
zi_spline=interp2(x,y,z,xi,yi,‘spline’);
zi_cubic=interp2(x,y,z,xi,yi,‘cubic’);
figure;
hold on;
subplot(321);
surf(x,y,z);
title(‘原始数据’);
subplot(322);
surf(xi,yi,zi_nearest);
title(‘临近点插值’);
subplot(323);
surf(xi,yi,zi_linear);
title(‘线性插值’);
subplot(324);
surf(xi,yi,zi_spline);
title(‘三次样条插值’);
subplot(325);
surf(xi,yi,zi_cubic);
title(‘三次多项式插值’);

三维插值