免费网站建设联系电话,乐清定制网站建设,科技太空讲座观后感,自己建网站流程要学什么文章目录 插值插值方法用Python解决插值问题 拟合最小二乘拟合数据拟合的Python实现 适用情况 处理由试验、测量得到的大量数据或一些过于复杂而不便于计算的函数表达式时#xff0c;构造一个简单函数作为要考察数据或复杂函数的近似 定义 给定一组数据#xff0c;需要确定满… 文章目录 插值插值方法用Python解决插值问题 拟合最小二乘拟合数据拟合的Python实现 适用情况 处理由试验、测量得到的大量数据或一些过于复杂而不便于计算的函数表达式时构造一个简单函数作为要考察数据或复杂函数的近似 定义 给定一组数据需要确定满足特定要求的曲线或曲面 插值如果所求曲线通过所给定有限个数据点 拟合要求所求曲线反映对象整体的变化态势得到简单实用的近似函数 插值
在一系列数据点对中一些数据点的函数值缺失因而希望通过已有数据点得到函数的近似表达式从而近似出确实数据点的函数值
从性质优良、便于计算的函数类{P(x)}中选择一个使得 P ( x i ) y i P(x_i) y_i P(xi)yi成立的P(x)作为f(x)的近似 x 0 , x 1 , . . . , x n x_0, x_1, ..., x_n x0,x1,...,xn成为插值节点 { P ( x ) } \{P(x)\} {P(x)}称为插值函数类 P ( x i ) y i ( i 0 , 1 , . . . , n ) P(x_i) y_i(i0, 1, ..., n) P(xi)yi(i0,1,...,n)称为插值条件 得到的 P ( x ) P(x) P(x)称为插值函数 f ( x ) f(x) f(x)称为被插值函数
一维插值方法一维Lagrange插值、分段线性插值、分段二次插值、牛顿插值和样条插值 二维插值方法二维样条插值
插值方法
Lagrange插值 P ( x ) ∑ i 0 n l i ( x ) y i P(x)\sum^n_{i0}l_i(x)y_i P(x)i0∑nli(x)yi 其中 l i ( x ) l_i(x) li(x)称为以 x 0 , x 1 , . . . , x n x_0, x_1, ..., x_n x0,x1,...,xn为节点的Lagrange插值基函数 l i ( x ) ∏ j 0 , j ≠ i n x − x j x i − x j l_i(x) \prod^n_{j0, j\neq i} \frac{x-x_j}{x_i-x_j} li(x)j0,ji∏nxi−xjx−xj
代码实现
def h(x,y,a):s 0.0for i in range(len(y)):t y[i]for j in range(len(y)):if i ! j:t * (a-x[j])(x[i]-x[j])s treturn s分段线性插值 用折线代替曲线 y f ( x ) y f(x) yf(x)其中 P ( x ) P(x) P(x)为 P ( x ) x − x i x i 1 − x i y i 1 x − x i 1 x i − x i 1 y i P(x) \frac{x-x_i}{x_{i1}-x_i}y_{i1} \frac{x-x_{i1}}{x_i-x_{i1}}y_i P(x)xi1−xix−xiyi1xi−xi1x−xi1yi 其中 x ∈ [ x i , x i 1 ] , i 0 , 1 , . . . , n − 1 x \in [x_i, x_{i1}], i0,1,..., n-1 x∈[xi,xi1],i0,1,...,n−1
分段二次插值 P ( x ) P(x) P(x)为一二次多项式即适用分段抛物线代替 y f ( x ) yf(x) yf(x)
牛顿插值 差分定义函数 f ( x ) f(x) f(x)等距节点 x i x 0 i h ( i 0 , 1 , . . . , n ) x_ix_0ih(i0, 1, ..., n) xix0ih(i0,1,...,n)一阶前向差分 Δ f i f i 1 − f i \Delta f_i f_{i1}-f_i Δfifi1−fi 高阶差分为差分的差分 Δ 0 f ( x ) f ( x ) \Delta^0 f(x) f(x) Δ0f(x)f(x) Δ m f ( x ) Δ m − 1 f ( x h ) − Δ m − 1 f ( x ) \Delta^m f(x) \Delta^{m-1}f(xh) - \Delta^{m-1}f(x) Δmf(x)Δm−1f(xh)−Δm−1f(x)
递归函数计算差分
def diff_forward(f, k, h, x):if k0: return f(x)else: return diff_forward(f, k-1, h, xh) - diff_forward(f, k-1, h, x)差商定义函数 f ( x ) f(x) f(x)相异节点 x 0 x 1 . . . x n x_0 x_1... x_n x0x1...xn 函数 f ( x ) f(x) f(x)关于节点 x i x_i xi x j x_j xj的一阶差商 f [ x i , x j ] f[x_i, x_j] f[xi,xj]有 f [ x i , x j ] f ( x i ) − f ( x j ) x i − x j f[x_i, x_j] \frac{f(x_i)-f(x_j)}{x_i-x_j} f[xi,xj]xi−xjf(xi)−f(xj) f ( x ) f(x) f(x)关于点 x i x_i xi x j x_j xj x k x_k xk的二阶差商有 f [ x i , x j , x k ] f [ x i , x j ] − f [ x j , x k ] x i − x k f[x_i, x_j, x_k] \frac{f[x_i, x_j]-f[x_j, x_k]}{x_i-x_k} f[xi,xj,xk]xi−xkf[xi,xj]−f[xj,xk] f ( x ) f(x) f(x)关于 x 0 , x 1 , . . . , x k x_0, x_1, ..., x_k x0,x1,...,xk的 k k k阶差商为 f [ x 0 , x 1 , . . . , x k ] f [ x 0 , x 1 , . . . , x k − 1 ] − f [ x 1 , x 2 , . . . , x k ] x 0 − x 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} f[x0,x1,...,xk]x0−xkf[x0,x1,...,xk−1]−f[x1,x2,...,xk]
代码示例计算 k 1 k1 k1个点对数据的 k k k阶差商
def diff_quo(xi[], fi[]):if len(xi)2 and len(fi)2:return (diff_quo(xi[:len(xi)-1],fi[:len(fi)-1])-diff_quo(xi[1:len(xi)],fi[1:len(fi)])) / float(xi[0]-xi[-1]) return (fi[0]- fi[1])/float(xi[0]-xi[1])Newton插值公式 一次Newton插值多项式 N 1 ( x ) f ( x 0 ) ( x − x 0 ) f [ x 0 , x 1 ] N_1(x)f(x_0)(x-x_0)f[x_0, x_1] N1(x)f(x0)(x−x0)f[x0,x1] 再根据各阶差商的定义可以得到 N n ( x ) N_n(x) Nn(x)即 f ( x ) f(x) f(x)的 n n n次插值多项式
样条插值 适用于对插值函数的光滑性有较高要求的问题 样条函数具有一定光滑性的分段多项式 给定 [ a , b ] [a,b] [a,b]的一个分划 Δ : a x 0 x 1 . . . x n b \Delta: ax_0 x_1 ... x_nb Δ:ax0x1...xnb 若 S ( x ) S(x) S(x)在各个小区间 [ x i , x i 1 ] ( i 0 , 1 , . . . , n − 1 ) [x_i, x_{i1}](i0, 1, ..., n-1) [xi,xi1](i0,1,...,n−1)上为 m m m次多项式且有 m − 1 m-1 m−1阶连续导数称 S ( x ) S(x) S(x)为关于分划 Δ \Delta Δ的 m m m次样条函数折线属于一次样条曲线
二维数据的双三次样条插值 考虑二维数据的插值时需要考虑到插值区域是否规则给定数据是有规律分布的还是散乱、随机分布的 当二维数据在规则区域上有规律分布时可以考虑用双三次样条插值即求解一个 S ( x ) S(x) S(x)满足对 x x x和 y y y都是三次的多项式
用Python解决插值问题
scipy.interpolatemodule有一维插值函数interp1d、二维插值函数interp2d和多维插值函数interpn
一维插值 interp1d(x, y, kindlinear) 说明kind参数取值为字符串指明插值方法取值包括linear线性插值、zero0阶样条插值、slinear1阶样条插值、quadratic2阶样条插值、cubic3阶样条插值 代码
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
xnp.arange(0,25,2)
ynp.array([12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13])
xnewnp.linspace(0, 24, 500) #插值点
f1interp1d(x, y);
y1f1(xnew);
f2interp1d(x, y,cubic);
y2f2(xnew)
plt.rc(font,size16);
plt.rc(font,familySimHei)
plt.subplot(121)
plt.plot(xnew, y1)
plt.xlabel(A分段线性插值)
plt.subplot(122)
plt.plot(xnew, y2)
plt.xlabel(B三次样条插值)
plt.savefig(figure7_4.png, dpi500)
plt.show()二维网格节点插值
思路原始数据为100x100网格节点的数据为提高精度适用双三次样条插值得到该区域上10x10网格节点的数据。把 0 ≤ x ≤ 1400 ∧ 0 ≤ y ≤ 1200 0 \leq x \leq 1400 \land 0 \leq y \leq 1200 0≤x≤1400∧0≤y≤1200 数据分为140x120个小矩形计算对应曲面面积每个矩形分为两个三角形再利用海伦公式求解 代码
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
import numpy as np
from numpy.linalg import norm
from scipy.interpolate import interp2d
znp.loadtxt(Pdata7_5.txt) #加载高程数据
xnp.arange(0,1500,100)
ynp.arange(1200,-100,-100)
finterp2d(x, y, z, cubic)
xnnp.linspace(0,1400,141)
ynnp.linspace(0,1200,121)
znf(xn, yn)
mlen(xn); nlen(yn); s0;
for i in np.arange(m-1):for j in np.arange(n-1):p1np.array([xn[i],yn[j],zn[j,i]])p2np.array([xn[i1],yn[j],zn[j,i1]])p3np.array([xn[i1],yn[j1],zn[j1,i1]])p4np.array([xn[i],yn[j1],zn[j1,i]])p12norm(p1-p2); p23norm(p3-p2); p13norm(p3-p1);p14norm(p4-p1); p34norm(p4-p3);L1(p12p23p13)/2;s1np.sqrt(L1*(L1-p12)*(L1-p23)*(L1-p13));L2(p13p14p34)/2; s2np.sqrt(L2*(L2-p13)*(L2-p14)*(L2-p34));sss1s2;
print(区域的面积为, s)
plt.rc(font,size16); plt.rc(text,usetexTrue)
plt.subplot(121); contrplt.contour(xn,yn,zn); plt.clabel(contr)
plt.xlabel($x$); plt.ylabel($y$,rotation90)
axplt.subplot(122,projection3d);
X,Ynp.meshgrid(xn,yn)
ax.plot_surface(X, Y, zn,cmapviridis)
ax.set_xlabel($x$); ax.set_ylabel($y$); ax.set_zlabel($z$)
plt.savefig(figure7_5.png,dpi500); plt.show()二维乱点插值 代码
from mpl_toolkits import mplot3d
import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata
xnp.array([129,140,103.5,88,185.5,195,105,157.5,107.5,77,81,162,162,117.5])
ynp.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])
xynp.vstack([x,y]).T
xnnp.linspace(x.min(), x.max(), 100)
ynnp.linspace(y.min(), y.max(), 100)
xng, yng np.meshgrid(xn,yn) #构造网格节点
zngriddata(xy, z, (xng, yng), methodnearest) #最近邻点插值
plt.rc(font,size16); plt.rc(text,usetexTrue)
axplt.subplot(121,projection3d);
ax.plot_surface(xng, yng, zn,cmapviridis)
ax.set_xlabel($x$); ax.set_ylabel($y$); ax.set_zlabel($z$)
plt.subplot(122); cplt.contour(xn,yn,zn,8); plt.clabel(c)
plt.savefig(figure7_6.png,dpi500); plt.show()拟合
最小二乘拟合
解决什么问题 已知一组二维数据即平面上 n n n个点 ( x i , y i ) ( i 1 , 2 , . . . , n ) (x_i, y_i)(i1, 2, ..., n) (xi,yi)(i1,2,...,n) x i x_i xi互不相同求函数 f ( x ) f(x) f(x)使得 f ( x ) f(x) f(x)在某种准则下与所有数据点最为接近即曲线拟合得最好 残差 δ i f ( x i ) − y i , i 1 , 2 , . . . , n \delta_if(x_i)-y_i, i1,2,...,n δif(xi)−yi,i1,2,...,n 最小二乘法使用的判定准则为残差的平和和最小即 a r g m i n J ∑ i 1 n ( f ( x i ) − y i ) 2 argmin \quad J\sum^n_{i1}(f(x_i)-y_i)^2 argminJi1∑n(f(xi)−yi)2 最终得到拟合函数 f ( x ) f ( x , a 1 , a 2 , . . . , a n ) f(x) f(x, a_1, a_2, ..., a_n) f(x)f(x,a1,a2,...,an) 根据参数 a 1 , a 2 , . . . , a n a_1, a_2, ..., a_n a1,a2,...,an线性与否最小二乘法分为线性最小二乘和非线性最小二乘
线性最小二乘法 给定线性无关的函数系 { ϕ k ( x ) ∣ k 1 , 2 , . . . , m } \{\phi_k(x)|k1,2,...,m\} {ϕk(x)∣k1,2,...,m} 若有拟合函数 f ( x ) ∑ k 1 m a k ϕ k ( x ) f(x) \sum^m_{k1}a_k \phi_k(x) f(x)∑k1makϕk(x)例如 f ( x ) a m x m − 1 a m − 1 x m − 2 . . . a 2 x a 1 f(x)a_mx^{m-1}a_{m-1}x^{m-2}...a_2xa_1 f(x)amxm−1am−1xm−2...a2xa1 或 f ( x ) ∑ k 1 m a k c o s ( k x ) f(x) \sum^m_{k1}a_k cos(kx) f(x)∑k1makcos(kx) 称 f ( x ) f ( x , a 1 , a 2 , . . . , a m ) f(x)f(x,a_1, a_2, ..., a_m) f(x)f(x,a1,a2,...,am)为关于参数 a 1 , a 2 , . . . , a m a_1,a_2,..., a_m a1,a2,...,am的线性函数 将 f ( x ) f(x) f(x)带入 J J J的计算根据 ∂ J ∂ a k 0 , k 1 , 2 , ⋯ , m \frac{\partial J}{\partial a_k}0,\quad k1,2,\cdots,m ∂ak∂J0,k1,2,⋯,m 即 ∑ i 1 n [ ( f ( x i ) − y i ) φ k ( x i ) ] 0 , k 1 , 2 , ⋯ , m \sum_{i1}^{n}\left[\left(f\left(x_{i}\right)-y_{i}\right) \varphi_{k}\left(x_{i}\right)\right]0, \quad k1,2, \cdots, m i1∑n[(f(xi)−yi)φk(xi)]0,k1,2,⋯,m 得到 形成一个关于 a 1 , a 2 , . . . , a m a_1,a_2,...,a_m a1,a2,...,am的线性方程组记号说明如下 则正规方程组改写为 R T R A R T Y R^TRAR^TY RTRARTY 当矩阵 R R R列满秩时 R T R R^TR RTR可逆此时正规方程组有唯一解即 A ( R T R ) − 1 R T Y A(R^TR)^{-1}R^TY A(RTR)−1RTY 非线性最小二乘拟合 当拟合函数不能以 ϕ k ( x ) \phi_k(x) ϕk(x)线性组合的形式构成时例如下列形式 使用同样的“最小化偏差”方法求解参数
拟合函数的选择 先做出数据的散点图从直观上判断应选用什么样的拟合函数 举个例子 若数据分布接近直线使用线性函数 f ( x ) a 1 x a 2 f(x)a_1xa_2 f(x)a1xa2 若数据分布接近抛物线使用二次多项式 f ( x ) a 1 x 2 a 2 x a 3 f(x)a_1x^2a_2xa_3 f(x)a1x2a2xa3 若数据分布是开始上升块随后逐渐变缓使用双曲线型函数或指数型函数即 若数据分布是开始下降较快随后逐渐变缓使用
数据拟合的Python实现
利用NumPy库中的多项式拟合函数polyfit或scipy.optimize模块中的curve_fit函数 polyfit的用法 代码展示
from numpy import polyfit, polyval, array, arange
from matplotlib.pyplot import plot,show,rc
x0arange(0, 1.1, 0.1)
y0array([-0.447, 1.978, 3.28, 6.16, 7.08, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2])
ppolyfit(x0, y0, 2) #拟合二次多项式
print(拟合二次多项式的从高次幂到低次幂系数分别为:,p)
yhatpolyval(p,[0.25, 0.35]); print(预测值分别为, yhat)
rc(font,size16)
plot(x0, y0, *, x0, polyval(p, x0), -); show()curve_fit的用法 调用格式 popt, pcov curve_fit(func, xdata, ydata) 参数说明func为拟合的函数xdata是自变量的观测值ydata是函数的观测值返回值popt是拟合的参数pcov是参数的协方差矩阵 代码展示
import numpy as np
from scipy.optimize import curve_fit
ylambda x, a, b, c: a*x**2b*xc
x0np.arange(0, 1.1, 0.1)
y0np.array([-0.447, 1.978, 3.28, 6.16, 7.08, 7.34, 7.66, 9.56, 9.48, 9.30, 11.2])
popt, pcovcurve_fit(y, x0, y0)
print(拟合的参数值为, popt)
print(预测值分别为, y(np.array([0.25, 0.35]), *popt))实例练习 代码
import numpy as np
from scipy.optimize import curve_fit
x0np.array([6, 2, 6, 7, 4, 2, 5, 9])
y0np.array([4, 9, 5, 3, 8, 5, 8, 2])
z0np.array([5, 2, 1, 9, 7, 4, 3, 3])
xy0np.vstack((x0, y0))
def Pfun(t, a, b, c):return a*np.exp(b*t[0])c*t[1]**2
popt, pcovcurve_fit(Pfun, xy0, z0)
print(abc的拟合值为, popt)代码
from mpl_toolkits import mplot3d
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
m200; n300
xnp.linspace(-6, 6, m); ynp.linspace(-8, 8, n);
x2, y2 np.meshgrid(x, y)
x3np.reshape(x2,(1,-1)); y3np.reshape(y2, (1,-1))
xynp.vstack((x3,y3))
def Pfun(t, m1, m2, s):return np.exp(-((t[0]-m1)**2(t[1]-m2)**2)/(2*s**2))
zPfun(xy, 1, 2, 3); zrz0.2*np.random.normal(sizez.shape) #噪声数据
popt, pcovcurve_fit(Pfun, xy, zr) #拟合参数
print(三个参数的拟合值分别为,popt)
znPfun(xy, *popt) #计算拟合函数的值
zn2np.reshape(zn, x2.shape)
plt.rc(font,size16)
axplt.axes(projection3d) #创建一个三维坐标轴对象
ax.plot_surface(x2, y2, zn2,cmapgist_rainbow)
plt.savefig(figure7_10.png, dpi500); plt.show()
