Mathematical Modeling - Interpolation and Fitting

1. Interpolation and fitting

1. What is interpolation and fitting?

 Interpolation and fitting is to construct a function through sample data. If the constructed function image is required to pass through all sample points, it is interpolation, otherwise it is fitting.

2. Commonly used one-dimensional interpolation methods

(1) piecewise linear interpolation

 Piecewise linear interpolation is to first divide the sample into n segments, and find a segmental polynomial P(x) to satisfy it on each segment: P i ( x ) = y i , i = 0 , 1... , n P_{i}(x)=y_{i} ,i=0,1...,n Pi​(x)=yi​,i=0,1...,n

(2) Spline interpolation

 Spline interpolation is to divide the sample into n segments, that is, to set (n+1) points. set point to x i , i = 0 , 1 , . . . , n x_{i},i=0,1,...,n xi​,i=0,1,...,n. have functions S ( x ) S(x) S(x) satisfies S i ( x ) S_{i}(x) Si​(x):
(1) in the interval [ x i , x i + 1 ] [x_{i},x_{i+1}] [xi​,xi+1​] is m m m degree polynomial
(2) in the interval [ x 0 , x n ] [x_{0},x_{n}] [x0​,xn​] on m − 1 m-1 m−1 order continuous derivative
then called S ( x ) S(x) S(x) is about this division m m spline function of degree m, m m The curve of the spline function of degree m is called a spline curve.
  It can be seen that the piecewise linear interpolation is actually a spline interpolation.
  One-dimensional interpolation in Python can be implemented with the interp1d function of the scipy.interpolate module.

3. Commonly used two-dimensional interpolation methods

(1) The interpolation node is grid data

 When the interpolation node is grid data, bicubic spline interpolation is generally used. In Python, it can be implemented with the interp2d function of the scipy.interpolate module.

This is the grid interpolation node

This is a non-grid interpolation node (if it is written in the form of the above table, there will only be z data on the diagonal of the table)

(2) The interpolation node is non-grid data

  Solved using the griddata function of the scipy.interpolate module.

4. Fitting

  Least squares fitting is a commonly used method for fitting curves in mathematics. There are fitting functions designed according to this principle in Python. For example, the curve_fit function of the scipy.optimize module uses nonlinear least squares fitting. and numpy's polyfit polynomial least squares fitting.
 Different from obtaining the interpolation function, the fitting function must be guessed before fitting. Guessing is best done on the basis of scatterplots of existing data. There are the following principles to choose the appropriate fitting function:

data trendsoptional function
tends to a straight line f ( x ) = a 1 x + a 2 f(x)=a_{1}x+a_{2} f(x)=a1​x+a2​
tends to parabola f ( x ) = a 1 x 2 + a 2 x + a 3 f(x)=a_{1}x^{2}+a_{2}x+a_{3} f(x)=a1​x2+a2​x+a3​
start rising fast then rise slowly f ( x ) = a 1 e − a 2 x f(x)=a_{1}e^{\frac{-a_{2}}{x}} f(x)=a1​ex−a2​or f ( x ) = x a 1 x + a 2 f(x)=\frac{x}{a_{1}x+a_{2}} f(x)=a1​x+a2​x​
start falling fast then fall slowly f ( x ) = a 1 e − a 2 x f(x)=a_{1}e^{-a_{2}x} f(x)=a1​e−a2​x or f ( x ) = 1 a 1 x + a 2 f(x)=\frac{1}{a_{1}x+a_{2}} f(x)=a1​x+a2​1​or f ( x ) = 1 a 1 x 2 + a 2 f(x)=\frac{1}{a_{1}x^{2}+a_{2}} f(x)=a1​x2+a2​1​

2. Interpolation function

1.interp1d function

1.grammar: scipy.interpolate.interp1d(x, y, kind='linear', axis=- 1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)
2.Function: interpolate one-dimensional function
3.Common parameters:
(1)x: real-valued vector representing the interpolation node
(2)y: N Dimensional real-valued array, the length of the interpolation direction is equal to x length
(3)kind: Interpolation method selection. linear(piecewise linear interpolation), nearest(nearest neighbor interpolation), slinear(single spline), quadratic(quadratic splines), cubic(cubic spline)
(4)axis: Specifies the interpolation direction, the default is the last axis. That is if y is a two-dimensional array, then the interpolation direction is horizontal ( axis==1)
4.Return value: returns the interpolation function f(x)

  Example
Get the interpolation function based on the following sample data and draw it

import numpy as np
from matplotlib import pyplot as plt
from scipy.interpolate import interp1d
# interpolation node
x=np.arange(0,25,2)
# The value corresponding to the interpolation node
y=np.array([12,9,9,10,18,24,28,27,25,20,18,15,13])
# interpolation point
x_new=np.linspace(0,24,500)
# Interpolation function, the first uses piecewise linear interpolation, the second uses cubic spline interpolation
f1=interp1d(x,y)
f2=interp1d(x,y,kind='cubic')
# draw
ax1=plt.subplot(1,2,1)
ax1.plot(x_new,f1(x_new))
ax2=plt.subplot(1,2,2)
ax2.plot(x_new,f2(x_new))

2.interp2d function

1.grammar: scipy.interpolate.interp2d(x, y, z, kind='linear', copy=True, bounds_error=False, fill_value=None)
2.Function: interpolate on a two-dimensional grid
3.Common parameters:
(1)x with y: Sample point coordinates x direction and y direction vector if x with y are multidimensional, then they are automatically flatten use later
(2)z: An array of interpolated point values, which will be compressed to one-dimensional before use
(3)kind: Interpolation method selection. The default is piecewise linear interpolation l(inear),Optional cubic splines ( cubic)or quintic spline interpolation ( quintic)
4.Return value: returns the interpolation function f(x,y)

  Example
Calculate the approximate value of the regional surface area based on the following elevation data, obtain the interpolation function and draw the regional contour line and three-dimensional surface map.

import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import norm #The role of the norm function is to find the matrix norm
from scipy.interpolate import interp2d
import warnings
from matplotlib import font_manager
my_font = font_manager.FontProperties(fname=r"C:\Windows\Fonts\Dengb.ttf")
warnings.filterwarnings('ignore')
# The value corresponding to the interpolation node
z=np.array([[1350,1370,1390,1400,1410,960,940,880,800,690,570,430,290,210,150],
          [1370,1390,1410,1430,1440,1140,1110,1050,950,820,690,540,380,300,210],
          [1380,1410,1430,1450,1470,1320,1280,1200,1080,940,780,620,460,370,350],
          [1420,1430,1450,1480,1500,1550,1510,1430,1300,1200,980,850,750,550,500],
          [1430,1450,1460,1500,1550,1600,1550,1600,1600,1600,1550,1500,1500,1550,1550],
          [950,1190,1370,1500,1200,1100,1550,1600,1550,1380,1070,900,1050,1150,1200],
          [910,1090,1270,1500,1200,1100,1350,1450,1200,1150,1010,880,1000,1050,1100],
          [880,1060,1230,1390,1500,1500,1400,900,1100,1060,950,870,900,936,950],
            [830,980,1180,1320,1450,1420,400,1300,700,900,850,810,380,780,750],
          [740,880,1080,1130,1250,1280,1230,1040,900,500,700,780,750,650,550],
          [650,760,880,970,1020,1050,1020,830,800,700,300,500,550,480,350],
          [510,620,730,800,850,870,850,780,720,650,500,200,300,350,320],
          [370,470,550,600,670,690,670,620,580,450,400,300,100,150,250]])
# interpolation node
x=np.arange(0,1500,100)
y=np.arange(1200,-100,-100)
#Use bicubic spline interpolation to get the interpolation function f
f=interp2d(x,y,z,kind='cubic')


# Set the subpoint xi=10i (i=0,1...140), yi=10j (j=0,1,...,120) the plane is divided into 140*120 small rectangles by this subpoint
xi=np.linspace(0,1400,141)
yj=np.linspace(1200,0,120)
# Interpolation of small rectangle vertices
zn=f(xi,yj)
# Calculation of area; the surface area corresponding to a small rectangle can be approximated as a quadrilateral in space. The area of ​​this quadrilateral can be regarded as the sum of the areas of the two triangles
# The triangle area is considered to be solved by Heron's formula: p=(a+b+c)/2 s=sqrt(p*(p-a)*(p-b)*(p-c))
s=0
for i in range(len(xi)-1):
    for j in range(len(yj)-1):
        p1=np.array([xi[i],yj[j],f(xi[i],yj[j])])
        p2=np.array([xi[i],yj[j+1],f(xi[i],yj[j+1])])
        p3=np.array([xi[i+1],yj[j],f(xi[i+1],yj[j])])
        p4=np.array([xi[i+1],yj[j+1],f(xi[i+1],yj[j])])
        a1,b1,c1=norm(p1-p2),norm(p1-p3),norm(p2-p3);p_1=(a1+b1+c1)/2;s1=np.sqrt(p_1*(p_1-a1)*(p_1-b1)*(p_1-c1))
        a2,b2,c2=norm(p4-p2),norm(p4-p3),norm(p3-p2);p_2=(a2+b2+c2)/2;s2=np.sqrt(p_2*(p_2-a2)*(p_2-b2)*(p_2-c2))
        s+=(s1+s2)
print('area{}cubic meter'.format(s[0]))


# Use plt.contour to draw contour line (contour line)
ax1=plt.subplot(1,2,1)
cont1=ax1.contour(xi,yj,zn)
# Use the plt.clabel function to add labels to contour lines
ax1.clabel(cont1)
plt.title('contour line')
# For a 3D surface plot, the subplot projection type is set to 3d. default 2d
ax2=plt.subplot(1,2,2,projection='3d')
X,Y=np.meshgrid(xi,yj)
ax2.plot_surface(X,Y,zn,cmap='viridis')
ax2.set_title('surface')
# full image title
plt.suptitle('Interpolate')


Regional area: 4880683.019677971 cubic meters

3.griddata function

1.grammar: scipy.interpolate.griddata(points, values, xi, method='linear', fill_value=nan, rescale=False)
2.Function: Do interpolation of non-grid interpolation node data
3.Common parameters:
(1)points: It is generally convenient to pass a tuple of data point coordinate vectors. Such as x is the interpolation node x a vector of values ​​for the direction, y is the interpolation node y vector of orientation values. So pass (x,y)
(2)values: is a vector of floating-point or complex numbers corresponding to the interpolation points
(3)xi: This parameter is the interpolation point. pass np.meshgrid function processing x,y The grid matrix is ​​obtained after interpolating the point vector in the direction and then passed as a tuple to xi More convenient
(4)kind: interpolation method. optional nearest,linear,cubic
4.Return value: return the interpolation point directly xi Corresponding interpolation array

  Example
Draw seabed topographic maps and contour maps based on water depth data

import  numpy as np
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
# Interpolate node vector data
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])
# Interpolation point x-direction and y-direction vector
xn=np.linspace(x.min(),x.max(),100)
yn=np.linspace(y.min(),y.max(),100)
# Construct grid nodes using vectors xn and yn
xng,yng=np.meshgrid(xn,yn)
# Use the griddata function to interpolate values ​​at grid nodes based on unstructured data
zn1=griddata((x,y),z,(xng,yng),method='nearest') # nearest neighbor interpolation
# draw 3D image
ax1=plt.subplot(131,projection='3d')
ax2=plt.subplot(133)
ax1.plot_surface(xng,yng,zn1,cmap='viridis')
contour=ax2.contour(xng,yng,zn1)
plt.clabel(contour)

3. Fitting function

1.polyfit function

1.grammar: numpy.polyfit(x, y, deg, rcond=None, full=False, w=None, cov=False)
2.Function: Least squares polynomial fitting, it seems that it can only do one-dimensional fitting
3.Common parameters:
(1)x: sample point x coordinate vector
(2)y: sample point y Coordinate vector or a two-dimensional array, each column represents the sample y coordinate vector, the columns share x coordinate vector
(3)deg: polynomial degree
4.Return Value: A vector of coefficients of the polynomial or a 2D array of the vectors of coefficients of the polynomial as columns

  Example
According to the sample data, a polynomial fitting curve is obtained and x = 0.25 with 0.35 x=0.25 and 0.35 Values ​​at x=0.25 and 0.35

import numpy as np
import matplotlib.pyplot as plt
# sample
x0=np.arange(0,1.1,0.1)
y0=np.array([-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.30,11.2])
# Use the polyfit function to fit the coefficients of the quadratic polynomial
p=np.polyfit(x0,y0,2)
# Use the polyval function to bring specified values ​​into polynomial calculations
yhat=np.polyval(p,np.array([0.25,0.35]))
print('x=0.25 predicted value of{:.2}  x=0.35 predicted value of{:.2}'.format(yhat[0],yhat[1]))
# plot the curve of the fitted function
plt.plot(x0,y0,'^')
plt.plot(x0,np.polyval(p,x0),'-')


2.curve_fit function

1.grammar: 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)
2.Function: Use nonlinear least square method to fit a function f
3.Common parameters:
(1)f: callable function f(x,...). He must treat the independent variables as f The first parameter of , the parameters in the function model are respectively as f the rest of the parameters. x can be viewed as a scalar or a vector, so we can define multidimensional functions
(2)xdata: sample point, if f is a one-dimensional function then xdata is a vector; if f yes k dimension function, x that is shape[0]==k array of . For example, my sample point is ( x,y,z),then my xdata can be written as np.vstack((x,y))
(3)ydata: A vector of function values ​​corresponding to the sample points
4.return value:
(1)popt: Optimal values for the parameters so that the sum of the squared residuals of f(xdata, *popt) - ydata is minimized.
(2)pcov: The estimated covariance of popt.

  Example
Fit function to sample points z = a e b x + c y 2 z=ae^{bx}+cy^{2} z=aebx+cy2

from scipy.optimize import curve_fit
import numpy as np
# Sample point data
x0=np.array([6,2,6,7,4,2,5,9])
y0=np.array([4,9,5,3,8,5,8,2])
z0=np.array([5,2,1,9,7,4,3,3])
# Fit function model 
def fun(var,a,b,c):
    return a*np.exp(b*var[0])+c*var[1]**2
# Curve fitting
popt,pcov=curve_fit(fun,np.vstack((x0,y0)),z0)
print("The parameters of the fitted function model are:{}".format(popt))

4. Drawing function

1.contour function

1.grammar: matplotlib.pyplot.contour(*args, data=None, **kwargs)
2.Function: draw contour line(contour line)
3.Common parameters:
(1)x: x vector of orientation values, satisfying len(x)==z.shape[0]
(2)y: y direction-worthy vector, satisfying len(y)==z.shape[1]
(3)z: The height value of the contour line, which is a two-dimensional array
(4)cmap:colormap can have cmap control, optional https://matplotlib.org/2.0.2/users/colormaps.html
(5)alpha: opacity, 0~1 between

2.plot_surface function

1.grammar: plot_surface(X, Y, Z, *, norm=None, vmin=None, vmax=None, lightsource=None, **kwargs)
2.Function: draw three-dimensional surface
3.Common parameters:
(1)X,Y,Z: A 2D array of sample points. pass np.meshgrid function can put x vector sum y Vectorize to 2D array X,Y
(2)cmap
(3)color

Tags: Python Mathematical Modeling numpy

Posted by russy on Fri, 30 Dec 2022 07:58:35 +0300