插值与拟合是数学建模中的一种基本的数据分析手段，被公认为建模中的常用算法之一。
插值
插值问题
例1.在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为：
12,9,9,10,18,24,28,27,25,20,18,15,13
推测中午一点温度，并作出24小时温度变化曲线图
例2.测得平板表面3*5网格点处的温度分别为：
82\quad81\quad80\quad82\quad84\\ 79\quad63\quad61\quad65\quad81\\ 84\quad84\quad82\quad85\quad86\\
作出平板表面的温度分布曲面z=f(x, y)的图形及等温线，并求出温度最高和最低点。
例3.在某海域测得一些点(x,y)处的水深z如下表，船的吃水深度为5英尺，在矩形区域(75,200)*(-50,150)内的哪些地方船要避免进入？
| x | y | z |
| 129 | 7.5 | 4 |
| 140 | 141.5 | 8 |
| 103.5 | 23 | 6 |
| 88 | 147 | 8 |
| 185.5 | 22.5 | 6 |
| 195 | 137.5 | 8 |
| 105 | 85.5 | 8 |
| 157.5 | -6.5 | 9 |
| 107.5 | -81 | 9 |
| 77 | 3 | 8 |
| 81 | 56.5 | 8 |
| 162 | -66.5 | 9 |
| 162 | 84 | 4 |
| 117.5 | -33.5 | 9 |
上述问题可归结为“已知函数在某区间(域)内若干点处的值,求函数在该区间(域)内其它点处的值”,这种问题适宜用插值方法解决。
一维插值问题可描述为：已知函数在x_0,x_1,...,x_n处的值y_0,y_1,...,y_n，求简单函数p(x)，使p(x_i)=y_i。
通常取p(x)为多项式，常用的插值方法有拉格朗日插值法和牛顿插值法。
拉格朗日插值法
拉格朗日插值公式（Lagrange Interpolation Formula）指的是在节点上给出节点基函数，然后做基函数的线性组合，组合系数为节点函数值的一种插值多项式。其定义如下：
对于某个多项式函数，已知有给定的k+1个取值点：
(x_0,y_0),...,(x_k,y_k)
其中x_j对应着自变量的位置，而y_j对应着函数在这个位置的取值。
假设任意两个不同的x_j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：
L(x)=\sum^k_{j=0}y_jl_j(x)其中每个l _{j}(x)为拉格朗日基本多项式（或称插值基函数），其表达式为：
l_j(x)=\prod_{i=0,i \not=j}^{k}\frac{x-x_i}{x_j-x_i}=\frac{x-x_0}{x_j-x_0}·...·\frac{x-x_{j-1}}{x_j-x_{j+1}}·\frac{x-x_{j+1}}{x_j-x_{j-1}}·...·\frac{x-x_k}{x_j-x_k}两式合并后为：
L(x)=\sum^k_{j=0}y_j\prod_{i=0,i \not=j}^{k}\frac{x-x_i}{x_j-x_i}拉格朗日基本多项式l_j(x)的特点是在x_j上取值为1，在其他的点x_i,i\not=j上取值为0。
牛顿插值法
牛顿多项式（英语：Newton Polynomial）是数值分析中一种用于插值的多项式，以英格兰数学家暨物理学家牛顿命名。其定义如下：
给定包含k+1个数据点的集合：
(x_0,y_0),...,(x_k,y_k)
如果对于∀i,j∈{0,…,k},i≠j，满足x_i \not= x_j，那么应用牛顿插值公式所得到的牛顿插值多项式为：
N(x)=\sum^k_{j=0}a_jn_j(x)其中每个n_j(x)为牛顿基本多项式（或称插值基函数），其表达式为：
n_j(x)=\prod_{i=0}^{j-1}(x-x_i),j>0,n_0(x)≡1系数a_j=[y_0,...,y_j]，而[y_0,...,y_j]表示差商。
差商表（高阶差商是两个低一阶差商的差商）
因此，牛顿多项式可以写作：
N(x)=[y_0]+[y_0,y_1](x-x_0)+...+[y_0,...,y_k](x-x_0)(x-x_1)...(x-x_{k-1})高次插值的Runge现象
在研究插值问题的初期，所有人都想当然地认为插值多项式的次数越高，插值精度越高。
Runge通过对一个例子的研究发现,上述结论仅仅在插值多项式的次数不超过七时成立;插值多项式的次数超过七时，插值多项式会出现严重的振荡现象，称之为Runge现象。
因此，在实际中不应使用七次以上的插值。
避免Runge现象的常用方法是:将插值区间分成若干小区间，在小区间内用低次(二次，三次)插值，即分段低次插值，如样条函数插值。
Python中的插值模块
插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。与拟合不同的是，要求曲线通过所有的已知数据。SciPy的interpolate模块提供了许多对数据进行插值运算的函数，范围涵盖简单的一维插值到复杂多维插值求解。当样本数据变化归因于一个独立的变量时，就使用一维插值；反之样本数据归因于多个独立变量时，使用多维插值。
计算插值有两种基本的方法：
1、对一个完整的数据集去拟合一个函数；
2、对数据集的不同部分拟合出不同的函数，而函数之间的曲线平滑对接。
第二种方法又叫做仿样内插法，当数据拟合函数形式非常复杂时，这是一种非常强大的工具。
当样本数据变化归因于一个独立的变量时，就使用一维插值；反之样本数据归因于多个独立变量时，使用多维插值。
一维插值的Python实现
一维插值interp1d类
class scipy.interpolate.interp1d(x, y, kind='linear', axis=-1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)这是一个类，用于完成一维数据的插值运算。
参数
| 参数 | 数据类型 | 意义 |
| x | array_like | 一维数据 |
| y | array_like | N维数据，其中插值维度的长度必须与x长度相同 |
| kind | str or int, optional | 给出插值的样条曲线的阶数 ‘zero’ 、’nearest’零阶 ‘slinear’ 、’linear’线性 ‘quadratic’ 、’cubic’二阶和三阶样条曲线，更高阶的曲线可以直接使用整数值指定 |
| axis | int, optional | 指定y中插值的轴，默认是y的最后一维 |
| copy | bool, optional | 如果True（default）类内置x和y的备份 |
| bounds_error | bool, optional | 如果True(Default)，在插值过程中超出x的范围就会报错ValueError; 如果False，超界的值由fill_value指定。 默认是True，除非fill_value=‘extrapolate’ |
| fill_value | array-like or (array-like, array_like) or “extrapolate”, optional | 如果是ndarray（或float），则此值将用于填充数据范围之外的请求点。如果未提供，则默认值为NaN。类阵列必须正确广播到非插值轴的尺寸。 如果是二元组，则第一个元素用作的填充值，第二个元素 用作的填充值。不是2元素元组的任何元素（例如list或ndarray，不考虑形状）都被视为一个类似数组的自变量，意味着将两个边界都用作 。x_newx[-1]below,above=fill_value,fill_value 如果“外推”，则数据范围外的点将被外推。 |
| assume_sorted | bool, optional | 如果为False，则x的值可以按任何顺序排列，并且将首先对其进行排序。如果为True，则x必须是单调递增值的数组。 |
属性
| 属性 | 意义 |
| fill_value | fill_value的值 |
方法
| 方法 | 意义 |
| _call_(x) | 评估插值逼近 |
代码示例(例1实现)
import numpy as np
from scipy import interpolate
import pylab as pl
x=np.linspace(0,24,13)
#x=[ 0. 2. 4. 6. 8. 10. 12. 14. 16. 18. 20. 22. 24.]
y=np.array([12,9,9,10,18,24,28,27,25,20,18,15,13])
pl.figure(figsize=(16,8),dpi=80)
xnew=np.linspace(0,24,1300) #重设精度
pl.plot(x,y,"ro",label='raw')
for kind in ["nearest","zero","slinear","quadratic","cubic"]:
#"nearest","zero"为阶梯插值；slinear 线性插值；"quadratic","cubic" 为2阶、3阶B样条曲线插值
f=interpolate.interp1d(x,y,kind=kind)
#f是一个函数，用这个函数就可以找插值点的函数值了
ynew=f(xnew)
pl.plot(xnew,ynew,label=str(kind))
pl.legend(loc="lower right")
pl.show()
单变量插值UnivariateSpline
介绍
interp1d不能外推运算（外插值），但是UnivariateSpline可以外插值。
调用方式如下：
class scipy.interpolate.UnivariateSpline(x,y,w=None,bbox=[None,None],k=3,s=None)- x,y是X-Y坐标数组
- w是每个数据点的权重值
- k为样条曲线的阶数
- s为平滑参数。
- s=0，样条曲线强制通过所有数据点
- s>0,满足\sum (w(y-spline(x)))^2 \leq s
代码示例
from scipy import interpolate
import numpy as np
import matplotlib.pyplot as plt
x1=np.linspace(0,10,20)
y1=np.sin(x1)
x2=np.linspace(0,20,200)
y2=np.sin(x2)+np.random.normal(loc=0,scale=1,size=len(x2))*0.2
sx1=np.linspace(0,12,100)
func1=interpolate.UnivariateSpline(x1,y1,s=0) #s=0,强制通过所有点
sy1=func1(sx1)
sx2=np.linspace(0,22,2000)
func2=interpolate.UnivariateSpline(x2,y2,s=8) #s>0,不强制通过所有点
sy2=func2(sx2)
pl.figure(figsize=(20,6),dpi=80)
plt.subplot(1,2,1)
plt.plot(x1,y1,'o')
plt.plot(sx1,sy1)
plt.subplot(1,2,2)
plt.plot(x2,y2,'o')
plt.plot(sx2,sy2)
plt.show()
网格点二维插值的Python实现
二维插值interp2d类
使用警告
介绍
class scipy.interpolate.interp2d(x, y, z, kind='linear', copy=True, bounds_error=False, fill_value=None)这是一个类，用于完成一维数据的插值运算。
参数
| 参数 | 数据类型 | 意义 |
| x,y | array_like | 定义数据点坐标的数组。数据点的坐标需要按照递增的顺序进行排序。 如果这些点位于一个规则的网格上，x可以指定列坐标，y指定行坐标，例如： x = [0,1,2]; y = [0,3]; z = [[1,2,3], [4,5,6]]否则，x和y必须指定每个点的完整坐标，例如： x = [0,1,2,0,1,2]; y = [0,0,0,3,3,3]; z = [1,4,2,5,3,6] |
| z | array_like | 要在数据点上插值的函数值。如果z是一个多维数组，在使用前假定Fortran排序(order=’F’)而将其压平。如果x和y指定了列和行的坐标，扁平化的z数组的长度是len(x)*len(y)；如果x和y指定了每个点的坐标，则len(z) == len(x) == len(y)。 |
| kind | {‘linear’, ‘cubic’, ‘quintic’}, optional | 要使用的样条插值的种类。默认使用‘linear’ |
| copy | bool, optional | 如果True（default）类内置x和y的备份 |
| bounds_error | bool, optional | 如果True(Default)，在插值过程中超出x的范围就会报错ValueError; 如果False，超界的值由fill_value指定。 默认是True，除非fill_value=‘extrapolate’ |
| fill_value | array-like or (array-like, array_like) or “extrapolate”, optional | 如果是ndarray（或float），则此值将用于填充数据范围之外的请求点。如果未提供，则默认值为NaN。类阵列必须正确广播到非插值轴的尺寸。 如果是二元组，则第一个元素用作的填充值，第二个元素 用作的填充值。不是2元素元组的任何元素（例如list或ndarray，不考虑形状）都被视为一个类似数组的自变量，意味着将两个边界都用作 。x_newx[-1]below,above=fill_value,fill_value 如果“外推”，则数据范围外的点将被外推。 |
注意
- interp2d()中，输入的x,y,z先用ravel()被转成了一维数组。
- func()的输入必须是一维的，输出是二维的。
- 插值的源数据必须是等距网格。不然的haul，运行不保存但结果不对。
- 以k作为样条次数，如果任何维度的点数少于k+1，将产生一个错误。
代码示例1(例2实现)
由于基础数据为5*3，如果选用的kind为cubic，即k=3，那k+1=4>3，将会产生报错，所以只能使用默认的linear方法
import numpy as np
from scipy import interpolate
import pylab as pl
import matplotlib as mpl
# X-Y轴分为5*3的网格
y,x=np.mgrid[0:2:3j,0:4:5j]
z=np.array([[82,81,80,82,84],
[79,63,61,65,81],
[84,84,82,85,86]])
#设置插值函数
func=interpolate.interp2d(x,y,z)
xi=np.linspace(0,4,500)
yi=np.linspace(0,2,300)
zi=func(xi,yi)
# 绘图
# 为了更明显地比较插值前后的区别，使用关键字参数interpolation='nearest'
# 关闭imshow()内置的插值运算。
pl.figure(figsize=(16,8),dpi=80)
pl.subplot(121)
im1=pl.imshow(np.flip(z,axis=0),extent=[0,4,0,2],cmap=mpl.cm.hot,interpolation='nearest',origin="lower")
# extent=[-1,1,-1,1]为x,y范围 favals为函数值
pl.colorbar(im1)
pl.subplot(122)
im2=pl.imshow(np.flip(zi,axis=0),extent=[0,4,0,2],cmap=mpl.cm.hot, interpolation='nearest', origin="lower")
pl.colorbar(im2)
pl.show()
代码示例2
这里由于基础数据更多，因此不再受k的限制
import numpy as np
from scipy import interpolate
import pylab as pl
import matplotlib as mpl
def func(x,y):
return (x+y)*np.exp(-5.0*(x**2+y**2))
pl.figure(figsize=(16,8),dpi=80)
# X-Y轴分为15*15的网格
y,x=np.mgrid[-1:1:15j,-1:1:15j]
fvals=func(x,y) # 计算每个网格点上的函数值 15*15的值
#三次样条二维插值
newfunc=interpolate.interp2d(x,y,fvals,kind='cubic')
# 计算100*100的网格上的插值
xnew=np.linspace(-1,1,100) #x
ynew=np.linspace(-1,1,100) #y
fnew=newfunc(xnew, ynew)
# 绘图
# 为了更明显地比较插值前后的区别，使用关键字参数interpolation='nearest'
# 关闭imshow()内置的插值运算。
pl.subplot(121)
im1=pl.imshow(fvals,extent=[-1,1,-1,1],cmap=mpl.cm.hot,interpolation='nearest',origin="lower")
# extent=[-1,1,-1,1]为x,y范围 favals为函数值
pl.colorbar(im1)
pl.subplot(122)
im2=pl.imshow(fnew,extent=[-1,1,-1,1],cmap=mpl.cm.hot, interpolation='nearest', origin="lower")
pl.colorbar(im2)
pl.show()
RectBivariateSpline类
介绍
以下的介绍直接摘自scipy的官方文档：
class scipy.interpolate.RectBivariateSpline(x, y, z, bbox=[None, None, None, None], kx=3, ky=3, s=0)Bivariate spline approximation over a rectangular mesh.Can be used for both smoothing and interpolating data.
矩形网格上的双变量样条曲线近似。可用于平滑和内插数据。
参数
x,y : array_like
1-D arrays of coordinates in strictly ascending order. Evaluated points outside the data range will be extrapolated.
严格按照升序排列的一维坐标数组。数据范围外的被评估点将被外推。
z : array_like
2-D array of data with shape (x.size,y.size).
形状为（x.size,y.size）的二维数据阵列。
bbox : array_like, optional
Sequence of length 4 specifying the boundary of the rectangular approximation domain, which means the start and end spline knots of each dimension are set by these values. By default, bbox=[min(x), max(x), min(y), max(y)].
长度为4的序列，指定矩形近似域的边界，这意味着每个维度的开始和结束花键结都由这些值来设置。默认情况下，bbox=[min(x), max(x), min(y), max(y)].
kx, ky : ints, optional
Degrees of the bivariate spline. Default is 3.
双变量样条曲线的度数。默认为3。
s : float, optional
Positive smoothing factor defined for estimation condition:
sum((z[i]-f(x[i],y[i]))**2,axis=0)<=s
where f is a spline function. Default is s=0, which is for interpolation.
为估计条件定义的正的平滑系数。其中f是一个样条函数。默认为s=0，用于插值。
注意
- 输入虽然标识为x,y，但事实上第一个参数是axis=0,第二个参数对应axis=1
- rectBivariateSpline中控制插值函数阶次的参数是kx,ky
- interp2d中控制插值方法的参数是kind ，kind=’cubic’ 表示三次插值
- 这两个函数，输入参数x，y均为x和y方向上的网格向量
- 在调用插值时，如果需要采样的点构成一个正立的矩形网格的话，那么输入x和y就输入向量就可以，它会输出 x*y 坐标空间上的2d插值结果。如果需要插值的是一系列离散的点，或者不构成正立网格，那么输入的x和y就是一系列点的坐标，它会输出(x,y)坐标空间上离散点的插值结果。指示插值点规则的参数是grid，默认是true的。
代码示例(例2实现)
import numpy as np
from scipy import interpolate
import pylab as pl
import matplotlib as mpl
# X-Y轴分为5*3的网格
x=np.linspace(0,4,5)
y=np.linspace(0,2,3)
z=np.array([[82,81,80,82,84],
[79,63,61,65,81],
[84,84,82,85,86]])
#设置插值函数
func=interpolate.RectBivariateSpline(x,y,z.T,kx=1,ky=1)
xi=np.linspace(0,4,500)
yi=np.linspace(0,2,300)
zi=func(xi,yi)
# 绘图
# 为了更明显地比较插值前后的区别，使用关键字参数interpolation='nearest'
# 关闭imshow()内置的插值运算。
pl.figure(figsize=(16,8),dpi=80)
pl.subplot(121)
im1=pl.imshow(np.flip(z,axis=0),extent=[0,4,0,2],cmap=mpl.cm.hot,interpolation='nearest',origin="lower")
# extent=[-1,1,-1,1]为x,y范围 favals为函数值
pl.colorbar(im1)
pl.subplot(122)
im2=pl.imshow(np.flip(zi.T,axis=0),extent=[0,4,0,2],cmap=mpl.cm.hot, interpolation='nearest', origin="lower")
pl.colorbar(im2)
pl.show()
RegularGridInterpolator类
使用警告
介绍
以下的介绍直接摘自scipy的官方文档：
class scipy.interpolate.RegularGridInterpolator(points, values, method='linear', bounds_error=True, fill_value=nan)Interpolation on a regular or rectilinear grid in arbitrary dimensions.
在任意尺寸的规则或直角网格上进行插值。
The data must be defined on a rectilinear grid; that is, a rectangular grid with even or uneven spacing. Linear, nearest-neighbor, spline interpolations are supported. After setting up the interpolator object, the interpolation method may be chosen at each evaluation.
数据必须定义在一个直角网格上；也就是说，一个具有均匀或不均匀间隔的矩形网格。支持线性、近邻、样条插值。在设置了插值器对象后，可以在每次评估时选择插值方法。
参数
points : tuple of ndarray of float, with shapes (m1, ), …, (mn, )
The points defining the regular grid in n dimensions. The points in each dimension (i.e. every elements of the points tuple) must be strictly ascending or descending.
定义n个维度的规则网格的点。每个维度中的点（即点元组中的每个元素）必须是严格意义上的上升或下降的。
values : array_like, shape (m1, …, mn, …)
The data on the regular grid in n dimensions. Complex data can be acceptable.
n维的规则网格上的数据。复杂的数据也可以接受。
method : str, optional
The method of interpolation to perform. Supported are “linear”, “nearest”, “slinear”, “cubic”, “quintic” and “pchip”. This parameter will become the default for the object’s __call__ method. Default is “linear”.
要执行的插值方法。支持 “linear”、”nearest”、”slinear”、”cubic”、”quintic “和 “pchip”。这个参数将成为对象的call方法的默认值。默认为 “linear”。
bounds_error : bool, optional
If True, when interpolated values are requested outside of the domain of the input data, a ValueError is raised. If False, then fill_value is used. Default is True.
如果是True，当插值要求在输入数据的领域之外时，会产生ValueError。如果是False，则使用fill_value。默认为True。
fill_value : float or None, optional
The value to use for points outside of the interpolation domain. If None, values outside the domain are extrapolated. Default is np.nan.
用于插值域外的点的值。如果没有，域外的值将被外推。默认为np.nan。
注意
1.Contrary to LinearNDInterpolator and NearestNDInterpolator, this class avoids expensive triangulation of the input data by taking advantage of the regular grid structure.
In other words, this class assumes that the data is defined on a rectilinear grid.
与LinearNDInterpolator和NearestNDInterpolator相反，该类通过利用常规网格结构避免了对输入数据进行昂贵的三角测量。
换句话说，这个类假定数据是在一个直角网格上定义的。
2.The ‘slinear’(k=1), ‘cubic’(k=3), and ‘quintic’(k=5) methods are tensor-product spline interpolators, where k is the spline degree, If any dimension has fewer points than k + 1, an error will be raised.
slinear'(k=1), ‘cubic'(k=3), and ‘quintic'(k=5)方法是样条插值，其中k是样条程度，如果任何维度的点数少于k+1，将产生一个错误。
3.If the input data is such that dimensions have incommensurate units and differ by many orders of magnitude, the interpolant may have numerical artifacts. Consider rescaling the data before interpolating.
如果输入的数据是这样的：尺寸的单位不相称，并且相差很多数量级，内插法可能会有数字伪影。考虑在内插之前重新调整数据的比例。
代码示例(例2实现)
有意思的是在速度上RegularGridInterpolator比interp2d慢了整整4倍，所以不是很明白为什么scipy的那帮家伙会选择直接扔掉interp2d（（
import numpy as np
from scipy import interpolate
import pylab as pl
import matplotlib as mpl
# X-Y轴分为5*3的网格
x=np.linspace(0,4,5)
y=np.linspace(0,2,3)
z=np.array([[82,81,80,82,84],
[79,63,61,65,81],
[84,84,82,85,86]])
#设置插值函数
func=interpolate.RegularGridInterpolator((x,y),z.T,method='slinear')
yi,xi=np.mgrid[0:2:300j,0:4:500j]
zi=func((xi,yi))
# 绘图
# 为了更明显地比较插值前后的区别，使用关键字参数interpolation='nearest'
# 关闭imshow()内置的插值运算。
pl.figure(figsize=(16,8),dpi=80)
pl.subplot(121)
im1=pl.imshow(np.flip(z,axis=0),extent=[0,4,0,2],cmap=mpl.cm.hot,interpolation='nearest',origin="lower")
# extent=[-1,1,-1,1]为x,y范围 favals为函数值
pl.colorbar(im1)
pl.subplot(122)
im2=pl.imshow(np.flip(zi,axis=0),extent=[0,4,0,2],cmap=mpl.cm.hot, interpolation='nearest', origin="lower")
pl.colorbar(im2)
pl.show()
网格点二维插值综合实战：地形图
问题提出
在某山区测得一些地点的高程如下表。平面区域为0<x<5600，0<y<4800。尝试画出该山区的地貌图和等高线图，并求出最高和最低点。
| 4800 | 1350 | 1370 | 1390 | 1400 | 1410 | 960 | 940 | 880 | 800 | 690 | 570 | 430 | 290 | 210 | 150 |
| 4400 | 1370 | 1390 | 1410 | 1430 | 1440 | 1140 | 1110 | 1050 | 950 | 820 | 690 | 540 | 380 | 300 | 210 |
| 4000 | 1380 | 1410 | 1430 | 1450 | 1470 | 1320 | 1280 | 1200 | 1080 | 940 | 780 | 620 | 460 | 370 | 350 |
| 3600 | 1420 | 1430 | 1450 | 1480 | 1500 | 1550 | 1510 | 1430 | 1300 | 1200 | 980 | 850 | 750 | 550 | 500 |
| 3200 | 1430 | 1450 | 1460 | 1500 | 1550 | 1600 | 1550 | 1600 | 1600 | 1600 | 1550 | 1500 | 1500 | 1550 | 1550 |
| 2800 | 950 | 1190 | 1370 | 1500 | 1200 | 1100 | 1550 | 1600 | 1550 | 1380 | 1070 | 900 | 1050 | 1150 | 1200 |
| 2400 | 910 | 1090 | 1270 | 1500 | 1200 | 1100 | 1350 | 1450 | 1200 | 1150 | 1010 | 880 | 1000 | 1050 | 1100 |
| 2000 | 880 | 1060 | 1230 | 1390 | 1500 | 1500 | 1400 | 900 | 1100 | 1060 | 950 | 870 | 900 | 936 | 950 |
| 1600 | 830 | 980 | 1180 | 1320 | 1450 | 1420 | 400 | 1300 | 700 | 900 | 850 | 810 | 380 | 780 | 750 |
| 1200 | 740 | 880 | 1080 | 1130 | 1250 | 1280 | 1230 | 1040 | 900 | 500 | 700 | 780 | 750 | 650 | 500 |
| 800 | 650 | 760 | 880 | 970 | 1020 | 1050 | 1020 | 830 | 800 | 700 | 300 | 500 | 550 | 480 | 350 |
| 400 | 510 | 620 | 730 | 800 | 850 | 870 | 850 | 780 | 720 | 650 | 500 | 200 | 300 | 350 | 320 |
| 0 | 370 | 470 | 550 | 600 | 670 | 690 | 670 | 620 | 580 | 450 | 400 | 300 | 100 | 150 | 250 |
| Y/X | 0 | 400 | 800 | 1200 | 1600 | 2000 | 2400 | 2800 | 3200 | 3600 | 4000 | 4400 | 4800 | 5200 | 5600 |
数据导入与初始化
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import matplotlib as mpl
from scipy import interpolate
import pylab as pl
import mpl_toolkits.mplot3d
# 导入数据
file='./data/网格插值_地形图.xlsx' # 文件目录
# 读取前处理：不设置列名行，行名为第一行，省略最后一行
data=pd.read_excel(file,header=None,index_col=0,skipfooter=1)
# 将数据转换为数组
heights=np.array(data) # 读入数据
heights=np.flip(heights,axis=0) # 考虑原始数据特点，进行垂直翻转构建原始数据点阵
# 建立横纵轴
x=np.linspace(0,5600,15) # 横轴坐标
y=np.linspace(0,4800,13) # 纵轴坐标
step=max(x)/(x.size-1) # 单位步长
# 3D作图数据预处理
area_3D=6000 #规定在6000*6000的平面上作图
# 这里需要对横纵轴进行重组，使得x=y，统一数值
heights_reshape=np.pad(heights,(0,(area_3D/step-heights.shape[1]+1).astype(np.int32)))
blank_row=np.zeros((heights_reshape.shape[1]-heights_reshape.shape[0],heights_reshape.shape[1]))
heights_reshape=np.r_[heights_reshape,blank_row]
y_reshape,x_reshape=np.mgrid[0:area_3D+step:step,0:area_3D+step:step]绘制原始数据草图
fig=plt.figure(figsize=(20,8),dpi=80)
# 图1：平面地形图
ax1=fig.add_subplot(121)
im1=ax1.imshow( heights, # 数据源
extent=[0,5600,0,4800], # 横纵轴的范围
cmap='terrain', # 使用分层设色地形图
origin="lower", # 图像原点设置在左下角
vmin=-500, # 图例最小值
vmax=2000 # 图例最大值
)
plt.colorbar(im1)
# 图2：立体地形图
ax2=fig.add_subplot(122,projection='3d')
im2=ax2.plot_surface(x_reshape, # 重组横轴
y_reshape, # 重组纵轴
heights_reshape, # 重组高度表
vmin=-500, # 图例最小值
vmax=2000, # 图例最大值
cmap='terrain' # 使用分层设色地形图
)
ax2.view_init(50,-80)
ax2.set_box_aspect((6,6,1))
plt.show()
插值
# 设置插值函数（两种方法均已提供）
func1=interpolate.RectBivariateSpline(x,y,heights.T)
func2=interpolate.RegularGridInterpolator((x,y),heights.T,method='cubic')
# 设置插值精度
rate=100 # 插值精度
xi=np.linspace(0,5600,14*rate+1) # x
yi=np.linspace(0,4800,12*rate+1) # y
step_i=max(xi)/(xi.size-1) # 单位步长
# 这里对横纵轴进行重组，使得xi=yi，统一数值
yi_reshape,xi_reshape=np.mgrid[0:area_3D+step_i:step_i,0:area_3D+step_i:step_i]
blank_row_i=np.zeros((xi.size-yi.size,xi.size))使用RectBivariateSpline进行插值
插值
heights_1=func1(xi,yi).T
# 高度重组
heights_1_reshape=np.pad(heights_1,(0,(area_3D/step_i-heights_1.shape[1]+1).astype(np.int32)))
blank_row=np.zeros((heights_1_reshape.shape[1]-heights_1_reshape.shape[0],heights_1_reshape.shape[1]))
heights_1_reshape=np.r_[heights_1_reshape,blank_row]绘制图像
fig=plt.figure(figsize=(18,8),dpi=80)
# 图1：平面地形图
ax1=fig.add_subplot(121)
im1=ax1.imshow( heights_1, # 数据源
extent=[0,5600,0,4800], # 横纵轴的范围
cmap='terrain', # 使用分层设色地形图
origin="lower", # 图像原点设置在左下角
vmin=-500, # 图例最小值
vmax=2000 # 图例最大值
)
plt.colorbar(im1)
# 图2：立体地形图
ax2=fig.add_subplot(122,projection='3d')
im2=ax2.plot_surface(xi_reshape, # 重组横轴
yi_reshape, # 重组纵轴
heights_1_reshape, # 重组高度表
vmin=-500, # 图例最小值
vmax=2000, # 图例最大值
cmap='terrain' # 使用分层设色地形图
)
ax2.view_init(50,-80)
ax2.set_box_aspect((6,6,1))
plt.show()
绘制等高线图
fig=plt.figure(figsize=(20,8),dpi=80)
levels=np.array([0,100,200,400,600,800,1000,1200,1600,1800])
np.pad(xi,(0,xi_reshape.shape[0]-xi.size)),np.pad(yi,(0,yi_reshape.shape[0]-yi.size))
# 图1：平面地形图
ax1=fig.add_subplot(121)
im1=ax1.contour(xi,yi,heights_1,levels=levels)
plt.colorbar(im1)
# 图2：立体地形图
ax2=fig.add_subplot(122)
im2=ax2.contourf(xi,
yi,
heights_1,
levels=levels,
vmin=-500, # 图例最小值
vmax=2000, # 图例最大值
cmap='terrain' # 使用分层设色地形图
)
plt.colorbar(im2)
plt.show()
使用RegularGridInterpolator进行插值
插值
yy,xx=np.mgrid[0:4800+step_i:step_i,0:5600+step_i:step_i]
heights_2=func2((xx,yy))
# 高度重组
heights_2_reshape=np.pad(heights_2,(0,(area_3D/step_i-heights_2.shape[1]+1).astype(np.int32)))
blank_row=np.zeros((heights_2_reshape.shape[1]-heights_2_reshape.shape[0],heights_2_reshape.shape[1]))
heights_2_reshape=np.r_[heights_2_reshape,blank_row]绘制图像
fig=plt.figure(figsize=(18,8),dpi=80)
# 图1：平面地形图
ax1=fig.add_subplot(121)
im1=ax1.imshow( heights_2, # 数据源
extent=[0,5600,0,4800], # 横纵轴的范围
cmap='terrain', # 使用分层设色地形图
origin="lower", # 图像原点设置在左下角
vmin=-500, # 图例最小值
vmax=2000 # 图例最大值
)
plt.colorbar(im1)
# 图2：立体地形图
ax2=fig.add_subplot(122,projection='3d')
im2=ax2.plot_surface(xi_reshape, # 重组横轴
yi_reshape, # 重组纵轴
heights_2_reshape, # 重组高度表
vmin=-500, # 图例最小值
vmax=2000, # 图例最大值
cmap='terrain' # 使用分层设色地形图
)
ax2.view_init(50,-80)
ax2.set_box_aspect((6,6,1))
plt.show()
绘制等高线图
fig=plt.figure(figsize=(20,8),dpi=80)
levels=np.array([0,100,200,400,600,800,1000,1200,1600,1800])
np.pad(xi,(0,xi_reshape.shape[0]-xi.size)),np.pad(yi,(0,yi_reshape.shape[0]-yi.size))
# 图1：平面地形图
ax1=fig.add_subplot(121)
im1=ax1.contour(xi,yi,heights_2,levels=levels)
plt.colorbar(im1)
# 图2：立体地形图
ax2=fig.add_subplot(122)
im2=ax2.contourf(xi,
yi,
heights_2,
levels=levels,
vmin=-500, # 图例最小值
vmax=2000, # 图例最大值
cmap='terrain' # 使用分层设色地形图
)
plt.colorbar(im2)
plt.show()
离散点二维插值的Python实现
bisplrep函数
介绍
以下介绍直接摘自Scipy官方文档：
scipy.interpolate.bisplrep(x, y, z, w=None, xb=None, xe=None, yb=None, ye=None, kx=3, ky=3, task=0, s=None, eps=1e-16, tx=None, ty=None, full_output=0, nxest=None, nyest=None, quiet=1)Find a bivariate B-spline representation of a surface.
Given a set of data points (x[i], y[i], z[i]) representing a surface z=f(x,y), compute a B-spline representation of the surface. Based on the routine SURFIT from FITPACK.
寻找一个曲面的双变量B-spline表示。
给出一组数据点（x[i], y[i], z[i]）代表一个曲面z=f(x,y), 计算该曲面的B-spline表示。基于FITPACK中的SURFIT程序。
参数
x, y, z : ndarray
Rank-1 arrays of data points.
数据点的Rank-1数组。
w : ndarray, optional
Rank-1 array of weights. By default w=np.ones(len(x)).
等级1的权重数组。默认情况下，w=np.ones(len(x))。
xb, xe : float, optional
End points of approximation interval in x. By default xb = x.min(), xe=x.max().
默认情况下，xb=x.min()，xe=x.max()，是近似区间的端点。
yb, ye : float, optional
End points of approximation interval in y. By default yb=y.min(), ye = y.max().
默认情况下，yb=y.min()，ye=y.max()，是近似区间的端点。
kx, ky : int, optional
The degrees of the spline (1 <= kx, ky <= 5). Third order (kx=ky=3) is recommended.
样条的阶数（1 <= kx, ky <= 5）。建议使用三阶（kx=ky=3）。
task : int, optional
If task=0, find knots in x and y and coefficients for a given smoothing factor, s. If task=1, find knots and coefficients for another value of the smoothing factor, s. bisplrep must have been previously called with task=0 or task=1. If task=-1, find coefficients for a given set of knots tx, ty.
如果task=0，找到x和y中的结和给定平滑系数s的系数。如果task=1，找到结和平滑系数s的另一个值的系数。如果任务=-1，为给定的结点tx, ty寻找系数。
s : float, optional
A non-negative smoothing factor. If weights correspond to the inverse of the standard-deviation of the errors in z, then a good s-value should be found in the range (m-sqrt(2*m),m+sqrt(2*m)) where m=len(x).
一个非负的平滑系数。如果权重对应于z中误差的标准差的倒数，那么一个好的s值应该在(m-sqrt(2m),m+sqrt(2m))范围内找到，其中m=len(x)。
eps : float, optional
A threshold for determining the effective rank of an over-determined linear system of equations (0 < eps < 1). eps is not likely to need changing.
确定过定线性方程组的有效等级的阈值（0 < eps < 1），eps不太可能需要改变。
tx, ty : ndarray, optional
Rank-1 arrays of the knots of the spline for task=-1
task=1的样条函数的Rank-1数组
full_output : int, optional
Non-zero to return optional outputs.
非零表示返回可选输出。
nxest, nyest : int, optional
Over-estimates of the total number of knots. If None then nxest = max(kx+sqrt(m/2),2*kx+3), nyest = max(ky+sqrt(m/2),2*ky+3).
高估了结的总数。如果值为None，那么nxest = max(kx+sqrt(m/2),2kx+3), nyest = max(ky+sqrt(m/2),2ky+3)。
quiet : int, optional
Non-zero to suppress printing of messages.
非零，以抑制信息的打印。
返回值
tck : array_like
A list [tx, ty, c, kx, ky] containing the knots (tx, ty) and coefficients (c) of the bivariate B-spline representation of the surface along with the degree of the spline.
一个列表[tx, ty, c, kx, ky]，包含曲面的双变量B-spline表示法的结（tx, ty）和系数（c），以及样条阶数。
fp : ndarray
The weighted sum of squared residuals of the spline approximation.
样条近似的加权残差平方和。
ier : int
An integer flag about splrep success. Success is indicated if ier<=0. If ier in [1,2,3] an error occurred but was not raised. Otherwise an error is raised.
一个关于splrep成功的整数标志。如果ier<=0，则表示成功。如果ier在[1,2,3]中，则发生错误但没有被提出。否则，一个错误被提出。
msg : str
A message corresponding to the integer flag, ier.
一个与整数标志ier相对应的信息。
注意
See bisplev to evaluate the value of the B-spline given its tck representation.
If the input data is such that input dimensions have incommensurate units and differ by many orders of magnitude, the interpolant may have numerical artifacts. Consider rescaling the data before interpolation.
参见bisplev来评估B-spline给定的tck表示法的值。
如果输入的数据是这样的：输入尺寸的单位不相称，并且相差很多数量级，那么插值可能有数字伪影。可以考虑在内插之前重新调整数据的比例。
bisplev函数
介绍
以下介绍直接摘自Scipy官方文档：
scipy.interpolate.bisplev(x, y, tck, dx=0, dy=0)Evaluate a bivariate B-spline and its derivatives.
Return a rank-2 array of spline function values (or spline derivative values) at points given by the cross-product of the rank-1 arrays x and y. In special cases, return an array or just a float if either x or y or both are floats. Based on BISPEV and PARDER from FITPACK.
评估一个双变量B-spline函数和它的导数。
在特殊情况下，如果x或y或两者都是浮点数，则返回一个数组或仅仅是一个浮点数的样条函数值（或样条导数值）。基于FITPACK的BISPEV和PARDER。
参数
x, y : ndarray
Rank-1 arrays specifying the domain over which to evaluate the spline or its derivative.
指定评估样条曲线或其导数的域的第1级数组。
tck : tuple
A sequence of length 5 returned by bisplrep containing the knot locations, the coefficients, and the degree of the spline: [tx, ty, c, kx, ky].
bisplrep返回的长度为5的序列，包含结点位置、系数和样条阶数：[tx, ty, c, kx, ky]。
dx, dy : int, optional
The orders of the partial derivatives in x and y respectively.
分别为x和y的偏导数的阶数。
返回值
vals : ndarray
The B-spline or its derivative evaluated over the set formed by the cross-product of x and y.
B-spline或其导数在x和y的交叉乘积形成的集合上进行评估。
注意
See bisplrep to generate the tck representation.
参见bisplrep来生成tck表示法。
使用bisplrep/bisplev实现的二维离散点插值
以下为例三的求解代码：
import numpy as np
from scipy import interpolate
import pylab as pl
import matplotlib as mpl
x=np.array([129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5])
y=np.array([7.5,141.5,23,147,22.5,137.5,85.5,-6.5,-81,3,56.5,-66.5,84,-33.5])
z=np.array([4,8,6,8,6,8,8,9,9,8,8,9,4,9])
xi=np.arange(75,201)
yi=np.arange(-50,151)
tck=interpolate.bisplrep(x,y,z,kx=1,ky=1,s=0)
zi=interpolate.bisplev(xi,yi,tck=tck)
fig=plt.figure(figsize=(18,8),dpi=80)
ax1=fig.add_subplot(121)
im1=ax1.imshow(zi,extent=[75,200,-50,150],cmap='Blues',origin="lower")
plt.colorbar(im1)
ax2=fig.add_subplot(122)
levels=np.array([-40,-20,-15,-10,-5,0,5,10,15,20,40])
im2=ax2.contourf(xi,yi,zi.T,levels=levels,cmap='RdBu')
plt.colorbar(im2)
plt.show()
拟合
拟合问题
例4.从1点到12点每隔1小时测量一次温度，测得的温度依次为：
5,8,9,15,25,29,31,30,22,25,27,24
(1)试估计每隔1/10小时的温度值，并做出温度变化图形。
(2)推测t=13.5时的温度。
例4的第1问是典型的插值问题，但第2问不宜用插值方法，因为13.5已超出所给数据范围，用插值函数外推插值区间外的数据会产生较大误差。
高次插值出现了Runge现象，外推t>12时的数据出现了巨大误差。
样条插值的第一个图表明样条插值可以很好地解决第一个问题，但第二个图显示，用样条插值外推t>12时的数据也会产生较大误差。
综上，第二问不宜用插值方法。
解决第2问的常用方法是，根据1到12点间的温度数据求出温度与时间之间的近似函数关系f(t)，由f(t)推断t =13.5时的温度。这种根据离散数据求数据间近似函数关系的问题称为曲线拟合问题。
拟合和插值的主要区别
(1)插值函数过已知点，而拟合函数不一定过已知点；
(2)插值主要用于求函数值，而拟合的主要目的是求函数关系,从而进行预测等进一步的分析。
当然，某些特定的问题既可以用插值也可以用拟合。
拟合的计算
曲线拟合需解决如下两个问题：
(1)线型的选择；
(2)线型中参数的计算。
线型的选择是拟合计算的关键和难点。通常主要根据专业知识和散点图确定线型。
线性拟合中参数的计算可采用最小二乘法，而非线性拟合参数的计算则要应用Gauss-Newton迭代法。
拟合的Python实现
numpy库多项式拟合
polyfit函数
官方文档地址（开摆了，懒得搬了）：numpy.polyfit — NumPy v1.24 Manual
numpy.polyfit(x,y,deg,rcond=None,full=False,w=None,cov=False)例4的代码实现
import numpy as np
import matplotlib.pyplot as plt
# 数据
time=np.array([1,2,3,4,5,6,7,8,9,10,11,12])
temp=np.array([5,8,9,15,25,29,31,30,22,25,27,24])
# 拟合
func=np.polyfit(time,temp,3) # 进行操作
poly=np.poly1d(func) # 拟合函数
print(poly)
print(f'f(13.5)={poly(13.5)}')
# 生成拟合值
x=np.arange(1,15,0.1)
y=poly(x)
# 作图
fig=plt.figure(figsize=(10,6),dpi=80)
ax1=fig.add_subplot(111)
im1=ax1.plot(time,temp,'o')
im2=ax1.plot(x,y)
plt.show()3 2 -0.006475 x - 0.3283 x + 7.128 x - 4.434 f(13.5)=16.034188034187995
scipy库自定义拟合
curve_fit函数
scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=True, bounds=(-inf, inf), method=None, jac=None, *, full_output=False, **kwargs)例4的代码实现
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
# 数据
time=np.array([1,2,3,4,5,6,7,8,9,10,11,12])
temp=np.array([5,8,9,15,25,29,31,30,22,25,27,24])
# 定义一个目标拟合函数，这里以三次函数为例
def func(x,a,b,c,d):
return a*(x**3)+b*(x**2)+c*x+d
# 拟合，popt返回为拟合目标系数，pcov为结果估计方差
popt,pcov=optimize.curve_fit(func,time,temp)
print(popt)
# 拟合数据
x=np.arange(1,15,0.1)
y=func(x,popt[0],popt[1],popt[2],popt[3])
# 作图
fig=plt.figure(figsize=(10,6),dpi=80)
ax1=fig.add_subplot(111)
im1=ax1.plot(time,temp,'o')
im2=ax1.plot(x,y)
plt.show()[-6.47500648e-03 -3.28282828e-01 7.12807563e+00 -4.43434343e+00]
常用拟合类CommonFitting
封装常用的拟合函数，形成一个自定义类，并提供了一些接口：
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
class CommonFitting:
# 指数逼近
def __Exponential1(self,x,a,b):
return(a*np.exp(b*x))
def __Exponential2(self,x,a,b,c,d):
return(a*np.exp(b*x)+c*np.exp(d*x))
# 傅里叶逼近
def __Fourier(self,x,a0,a1,b1,w):
return(a0+a1*np.cos(x*w)+b1*np.sin(x*w))
# 高斯逼近
def __Gaussian(self,x,a1,b1,c1):
return(a1*np.exp(-((x-b1)/c1)**2))
# 多项式逼近
def __Linear(self,x,a,b):
return(a*x+b)
def __Cubic(self,x,a,b,c,d):
return(a*x**3+b*x**2+c*x+d)
def __Quadratic(self,x,a,b,c,d,e,f):
return(a*x**5+b*x**4+c*x**3+d*x**2+e*x+f)
# 幂逼近
def __Pow1(self,x,a,b):
return(a*x**b)
def __Pow2(self,x,a,b,c):
return(a*x**b+c)
# 正弦曲线逼近
def __SinFunc(self,x,a1,b1,c1):
return(a1*np.sin(b1*x+c1))
# 韦伯逼近
def __Weibull(self,x,a,b):
return(a*b*x**(b-1)*np.exp(-a*x**b))
# 拟合数据
def __doFitting(self,x):
if self.method=='Exponential1':
return self.__Exponential1(x,self.popt[0],self.popt[1])
elif self.method=='Exponential2':
return self.__Exponential2(x,self.popt[0],self.popt[1],self.popt[2],self.popt[3])
elif self.method=='Fourier':
return self.__Fourier(x,self.popt[0],self.popt[1],self.popt[2],self.popt[3])
elif self.method=='Gaussian':
return self.__Gaussian(x,self.popt[0],self.popt[1],self.popt[2])
elif self.method=='Linear':
return self.__Linear(x,self.popt[0],self.popt[1])
elif self.method=='Cubic':
return self.__Cubic(x,self.popt[0],self.popt[1],self.popt[2],self.popt[3])
elif self.method=='Quadratic':
return self.__Quadratic(x,self.popt[0],self.popt[1],self.popt[2],self.popt[3],self.popt[4],self.popt[5])
elif self.method=='Pow1':
return self.__Pow1(x,self.popt[0],self.popt[1])
elif self.method=='Pow2':
return self.__Pow2(x,self.popt[0],self.popt[1],self.popt[2])
elif self.method=='SinFunc':
return self.__SinFunc(x,self.popt[0],self.popt[1],self.popt[2])
elif self.method=='Pow1':
return self.__Weibull(x,self.popt[0],self.popt[1])
def __init__(self,x_data,y_data,method,target):
# 参数解释:x_data/y_data-原始数据;method-选择的拟合函数;target-目标拟合范围
# 导入数据
self.x=x_data # 横轴数据值
self.y=y_data # 纵轴数据值
self.xi=target # 目标x轴
# 逼近方法字典
self.method=method
methods={'Exponential1':self.__Exponential1, # a*np.exp(b*x)
'Exponential2':self.__Exponential2, # a*np.exp(b*x)+c*np.exp(d*x)
'Fourier':self.__Fourier, # a0+a1*np.cos(x*w)+b1*np.sin(x*w)
'Gaussian':self.__Gaussian, # a1*np.exp(-((x-b1)/c1)**2)
'Linear':self.__Linear, # a*x+b
'Cubic':self.__Cubic, # a*x**3+b*x**2+c*x+d
'Quadratic':self.__Quadratic, # a*x**5+b*x**4+c*x**3+d*x**2+e*x+f
'Pow1':self.__Pow1, # a*x**b
'Pow2':self.__Pow2, # a*x**b+c
'SinFunc':self.__SinFunc, # a1*np.sin(b1*x+c1)
'Weibull':self.__Weibull, # a*b*x**(b-1)*np.exp(-a*x**b)
}
# 选择拟合方式
if type(self.method)==str: # 进行自定义函数拟合
self.func=methods[self.method] # 选用的拟合方法
# 拟合目标系数,目标相对方差
self.popt,self.pcov=optimize.curve_fit(self.func,self.x,self.y)
self.yi=self.__doFitting(x=self.xi) # 目标y轴
elif type(self.method)==int: # 进行多项式拟合
self.popt=np.polyfit(self.x,self.y,self.method) # 拟合函数
self.func=np.poly1d(self.popt) # 拟合多项式
self.yi=self.func(self.xi) # 目标y轴
else:
raise TypeError("The type of method only can be 'str' or 'int'")
def __str__(self):
def addPlus(x):
if x>=0:
return(f'+{round(x,4)}')
else:
return(f'{round(x,4)}')
# 这里直接作打印拟合函数的操作
if type(self.method)==int:
res=''
for i in range(0,len(self.popt)):
if len(self.popt)-i-1==1:
res+=f'{addPlus(self.popt[i])}x'
elif len(self.popt)-i-1==0:
res+=f'{addPlus(self.popt[i])}'
else:
res+=f'{addPlus(self.popt[i])}x^{len(self.popt)-i-1}'
return res
else:
raise SyntaxError("Print custom functions is not completed,please try 'print(a.popt)'!")
# 搜寻拟合值
def find(self,x):
if type(self.method)==str:
return self.__doFitting(x)
else:
return self.func(x)
# 展示拟合图像
def show(self,size=(16,8)):
fig=plt.figure(figsize=size,dpi=80)
ax1=fig.add_subplot(111)
im1=ax1.plot(self.x,self.y,'o')
im2=ax1.plot(self.xi,self.yi)
plt.show()