数学建模算法与应用|插值方法(一维插值、三次样条插值、二维插值的matlab自带函数，python实现/作图) matlab|python

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

一维插值函数

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

‘nearest’	最近项插值
‘linear’	线性插值
‘spline’	立方样条插值
‘cubic’	立方插值

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

三次样条插值

y = interp1(x0, y0, x, ‘spline’)
y = spline(x0, y0, x)
pp = csape(x0, y0, conds)
pp = csape(x0, y0, conds, valconds); y = fnal(pp, x)

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

‘complete’	一阶导数
‘not-a-knot’	非扭结条件（没有边界条件）
‘periodic’	周期条件
'second‘	二阶导数

valconds默认导数值 [ 0, 0 ]

例1

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

x	0	3	5	7	9	11	12	13	14	15
y	0	1.2	1.7	2.0	2.1	2.0	1.8	1.2	1.0	1.6

matlab求插值

x0=[035791112131415]; y0=[01.21.72.02.12.01.81.21.01.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')

数学建模算法与应用|插值方法(一维插值、三次样条插值、二维插值的matlab自带函数，python实现/作图)

文章图片

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<...
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=[636697624478450 698712630478420 680674598412400 662626552334310]; 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')

"最高点的地址：""170180" "最高点的高程：""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*4f = 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自带函数，python实现/作图)】

文章图片

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);