插值方法(一维插值、三次样条插值、二维插值的matlab自带函数，python实现/作图)

• 数模比赛中，常常需要根据已知的函数点进行数据、模型的处理和分析，而有时候现有的数据是极少的，不足以支撑分析的进行，这时就需要使用一些数学的方法，“模拟产生”一些新的单又比较靠谱的值来满足需求，这就是插值的作用。
• 插值法在数值分析课程中有详细介绍。

一维插值函数

• y = interp1(x0, y0, x, ‘menthod’)
• **method **指定插值的方法，默认为线性插值。其值可为：

• 当 x0 为等距时可以用快速插值法，使用快速插值法的格式为‘*nearest’, ‘*linear’, ‘*spline’, ‘*cubic’

三次样条插值

1. y = interp1(x0, y0, x, ‘spline’)
2. y = spline(x0, y0, x)
3. pp = csape(x0, y0, conds)
4. pp = csape(x0, y0, conds, valconds); y = fnal(pp, x)

• 对于三次样条插值，提倡使用函数 csape 。
• condas: csape默认边界条件为Lagrange边界条件，其值可为：

• valconds默认导数值 [ 0, 0 ]

例1

• 给出下表数据，需要得到 x 坐标每改变0.1时的y坐标，并画出曲线。用分段线性和三次样条插值方法计算。

matlab求插值

x0=[0   3   5   7   9   11   12   13   14  15];
y0=[0  1.2  1.7  2.0  2.1  2.0  1.8  1.2   1.0  1.6];
x = 0:0.1:15;
y1 = interp1(x0,y0,x);%默认线性插值
y2 = interp1(x0,y0,x,'spline');%三次样条插值
pp1=csape(x0,y0);%边界条件为Lagrange条件
y3=fnval(pp1,x);

subplot(1,3,1)
plot(x0,y0,'+',x,y1)
title('Piecewise linear')

subplot(1,3,2)
plot(x0,y0,'+',x,y2)
title('Spline1')

subplot(1,3,3)
plot(x0,y0,'+',x,y3)
title('Spline2')

python求插值

import numpy as np
from scipy import interpolate #插值
import matplotlib.pyplot as plt #Pyplot 是 Matplotlib 的子库，提供了和 MATLAB 类似的绘图 API。Pyplot 是常用的绘图模块，能很方便让用户绘制 2D 图表。
x0 = np.array([0, 3, 5, 7, 9, 11, 12, 13, 14, 15])
y0 = np.array([0, 1.2, 1.7, 2.0, 2.1, 2.0, 1.8, 1.2, 1.0, 1.6])
x1 = np.arange(0, 15, 0.1)#numpy.arange(start, stop, step, dtype)  np.linspace(start, stop, num=50, endpoint=True, retstep=False, dtype=None)
plt.figure(figsize=(8, 6))#创建一个8*6点（point）的点图

###将图作在一张图上
for method in [ "slinear", "cubic"]:  # 插值方式
    f = interpolate.interp1d(x0, y0, kind=method) #一维数据的插值运算可以通过方法 interp1d() 完成。kind=method指插值方法
    y1 = f(x1)#得到插值结果
    plt.plot(x1, y1, label=method)
plt.plot(x0,y0,'o',label="datas")
plt.title('Interpolation')
plt.legend(loc="lower right")#图例的位置
plt.show()

###3张子图
i=1
for method in [ "slinear", "cubic"]:  # 插值方式
    f = interpolate.interp1d(x0, y0, kind=method) #一维数据的插值运算可以通过方法 interp1d() 完成。kind=method指插值方法
    y1 = f(x1)#得到插值结果
    plt.subplot(1, 3, i)
    plt.plot(x1, y1)
    plt.title(method)
    i=i+1
plt.subplot(1, 3, i)
plt.plot(x0,y0,'o-')
plt.title("Piecewise linear")
plt.suptitle('Interpolation')
plt.show()

二维插值

插值节点为网格节点

• 已知 m × n 个 节 点 ： ( x i , y j , z i j ) ( i = 1 , . . . , m , j = 1 , . . . , n ) ， 且 x 1 < . . . < x m , y 1 < . . . < y n m\times n个节点：(x_i,y_j,z_{ij})(i=1,...,m,j=1,...,n)，且x_1<...<x_m,y_1<...<y_n m×n个节点：(xi​,yj​,zij​)(i=1,...,m,j=1,...,n)，且x1​<...<xm​,y1​<...<yn​。求点 ( x , y ) 处 的 插 值 。 (x,y)处的插值。 (x,y)处的插值。
• z = interp2( x0, y0, z0, x, y, ’ method ') 注意：x0, y0分别是m维和n维向量，z0是n×m矩阵。x是M维行向量，y是N维列向量，z是N×M矩阵
• 三次样条插值** pp = csape( {x0, y0}, z0, conds, valconds), z=fnval(pp, {x, y})** 注意：x0, y0分别是m维和n维向量，z0是m×n矩阵。x是M维向量，y是N维向量，z是M×N矩阵

例2

• 在一丘陵地带测量高程， x x x和 y y y方向每隔100m 测一个点，得到高程表如下，试插值一个曲面，确定合适的模型，并由此找到最高点和该点的高程。

matlab求解

clear,clc
x=100:100:500;
y=100:100:400;
z=[636    697    624    478   450  
   698    712    630    478   420
   680    674    598    412   400
   662    626    552    334   310];

xi=100:10:500;yi=100:10:400;
pp=csape({x,y},z');
cz=fnval(pp,{xi,yi});
[i,j]=find(cz==max(max(cz)));
disp(["最高点的地址：" num2str([xi(i),yi(j)])]);
disp(["最高点的高程：" num2str(cz(i,j))]);
[X,Y]=meshgrid(xi,yi);
surf(X,Y,cz')

"最高点的地址："    "170  180"
"最高点的高程："    "720.6252"

python求解

import numpy as np
from scipy import interpolate #插值
import matplotlib.cm as cm #曲面的配色
import matplotlib.pyplot as plt #绘图
from mpl_toolkits.mplot3d import Axes3D #3D图

x0, y0 = np.mgrid[100:500:5j, 100:400:4j]#生成网格图5*4,5j代表5段，没有j代表间距为5
z0 = np.array([636, 697, 624, 478, 450,
               698, 712, 630, 478, 420,
               680, 674, 598, 412, 400,
               662, 626, 552, 334, 310]).reshape((4, 5))
z0=z0.T#注意z是5*4

f = interpolate.interp2d(x0, y0, z0, kind='cubic')#由样本点生成三次样条插值
x1 = np.linspace(100, 500, 100)#np.linspace(start, stop, num=50)
y1 = np.linspace(100, 400, 100)
z1 = f(x1, y1)#插值结果

x1, y1 = np.meshgrid(x1, y1)
fig = plt.figure()
ax=Axes3D(fig)#3D对象
#ax=fig.add_subplot(111,projection='3d')
#ax= plt.subplot(1, 2, 1, projection='3d')
surf = ax.plot_surface(x1, y1, z1, rstride=1, cstride=1,
                         cmap=cm.coolwarm, linewidth=0.5, antialiased=True)#rstride行之间的跨度，cstride列之间的跨度cmap颜色映射表，默认是“rainbow”
plt.title('Interpolation-2D(The peak: {:3f})'.format(np.max(z1)))
ax.scatter(x0,y0,z0,c='r')#三维散点图
ax.contour(x1, y1, z1,zdir='z',offset=300)#投影到z平面的z=300上
ax.set_zlabel('Z')  # 坐标轴
ax.set_ylabel('Y')
ax.set_xlabel('X')
plt.colorbar(surf, shrink=0.5, aspect=5)
plt.show()

#最高点处的地址
x1[np.where(z1==np.max(z1))] #array([144.44444444])
y1[np.where(z1==np.max(z1))] #array([200.])

插值节点为散乱节点

• 已知 n n n个节点， ( x i , y i , z i ) , i = 1 , 2 , . . . , n (x_i,y_i,z_i),i=1,2,...,n (xi​,yi​,zi​),i=1,2,...,n，求点 ( x , y ) (x,y) (x,y)处的插值 z z z。
• ZI = griddata(x, y, z, XI, YI ) 注意：x y z 均为n维向量，指明所给数据点的横坐标、纵坐标和竖坐标；向量XI YI是给定的网格点的横坐标和纵坐标，一个是行向量，另一个是列向量；返回值ZI 为网格（XI,YI）处的函数值。

例3

• 下表是海底水深数据，在适当矩形区域内画出海底曲面的图形。

matlab求解

clc, clear
clc, clear
x=[129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5];
y=[7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5];
z=-[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
xmm=minmax(x);  %求x的最小值和最大值
ymm=minmax(y);  %求y的最小值和最大值
xi=xmm(1):xmm(2);
yi=ymm(1):ymm(2);
zi1=griddata(x,y,z,xi,yi','cubic'); %立方插值
zi2=griddata(x,y,z,xi,yi','nearest'); %最近点插值
zi=zi1;  %立方插值和最近点插值的混合插值的初始值
zi(isnan(zi1))=zi2(isnan(zi1));  %把立方插值中的不确定值换成最近点插值的结果
subplot(1,2,1), plot(x,y,'*');
subplot(1,2,2), mesh(xi,yi,zi);