数学建模算法与应用 插值与拟合

导入包

import numpy as np

import matplotlib as mpl
import matplotlib.pyplot as plt
#中文显示问题
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False

from scipy.interpolate import CubicSpline
import scipy.integrate as si

习题

习题5.1

  在区间$[0,10]$上等间距取1000个点$x_i(i=1,2,\cdots ,1000)$，计算在这些点$x_i$处函数$g(x)=\frac{(3x^2+4x+6)sinx}{x^2+8x+6}$的函数值$y_i$，利用观测点$(x_i,y_i)(i=1,2,\cdots ,1000)$，求三次样条插值函数$\hat{g}(x)$，画出插值函数$\hat{g}(x)$的图形，并求积分${\int }_0^{10}g(x)dx$和${\int }_0^{10}\hat{g}(x)dx$

算法设计

# 定义g(x)
def gx(x):
    return (3*x**2+4*x+6)*np.sin(x)/(x**2+8*x+6)
x = np.linspace(0,10,1000)
y = gx(x)
# 求三次样条插值函数
cs = CubicSpline(x, y)
# 画出插值函数图像
plt.figure(dpi = 600)
plt.plot(x,cs(x))
# 求插值函数定积分
result = si.quad(cs, 0, 10)  # 函数 起点 终点
print(result)
# 求原函数定积分
result = si.quad(gx, 0, 10)  # 函数 起点 终点
print(result)

输出：

(2.243022802688952, 3.0531953334604206e-11)
(2.243022802741294, 2.0523107847822298e-11)

习题5.3

  已知当温度为$T=[700,720,740,760,780]$时，过热蒸汽体积的变化为$V=[0.0977,0.1218,0.1406,0.1551,0.1664]$，分别采用线性插值和三次样条插值求解T = 750,770时的体积变化，并在一个图形界面中画出线性插值函数和三次样条插值函数的图形。

算法设计

t = np.arange(700,800,20)
v = np.array([0.0977,0.1218,0.1406,0.1551,0.1664])
# 线性插值
f = inter.interp1d(t,v,kind ="linear")
# 求结果
v1 = f([750,770])
print('线性插值：',v1)
# 三次样条插值
cs = CubicSpline(t, v)
# 求结果
v2 = cs([750,770])
print('三次样条插值：',v2)
# 绘制图形
plt.figure(dpi = 600,figsize = (7,4))
plt.plot(t,v,'*-')
plt.plot(t,cs(t))
plt.legend(['线性插值','三次样条插值'])
plt.savefig('函数图像.png')

输出：

线性插值： [0.14785 0.16075]
三次样条插值： [0.14832031 0.16107969]

习题5.4

  某种合金的含铅量百分比为$p($，其溶解温度为$\theta {(}^{\circ }C)$，由试验测得$p$与$\theta$之间的经验公式$\theta =ap+b$。

p	36.9	46.7	63.7	77.8	84.0	87.5
$\theta$	181	197	235	270	283	292

算法设计

from sklearn.linear_model import LinearRegression
# 自变量
p = np.array([36.9,46.7,63.7,77.8,84.0,87.5])
# 因变量
y = np.array([181,197,235,270,283,292])
lin_reg = LinearRegression()#构建线性回归模型对象
# 使用np.newaxis增加一个维度
lin_reg.fit(p[:,np.newaxis],y[:,np.newaxis])
# 输出k值与b值（y = kx+b）
lin_reg.coef_,lin_reg.intercept_

输出：

(array([[2.23370015]]), array([95.35241998]))

习题5.5

   多项式$f(x)=a_3{x}^3+a_2{x}^2+a_{1}x+a_0$，取$a_3=8，a_2=5，a_1=2，a_0=-1$，在$[-6,6]$上等步长取$100$个点作为$x$的观测值，计算对应的函数值作为$y$的观测值；把得到的观测值记作$(x_i,y_i)，i=1,2,\cdots ,100$
  （1）利用观测值$(x_i,y_i)，i=1,2,\cdots ,100$，拟合三次多项式
  （2）把每个$y_i$加上白噪点，即加上一个服从标准正态分布的随机数，把得到的数据记作${\tilde{y}}_i$，利用$(x_i,{\tilde{y}}_i)，i=1,2,\cdots ,100$，拟合三次多项式

算法设计1

from sklearn.preprocessing import PolynomialFeatures
def fx(x):
    return 8 * x ** 3 + 5 * x ** 2 + 2 * x - 1
x = np.linspace(-6,6,100).reshape(-1, 1)
# 升维
x_0 = PolynomialFeatures(degree=3).fit_transform(x)
lin_reg = LinearRegression()#构建线性回归模型对象
lin_reg.fit(x_0,fx(x))
lin_reg.coef_,lin_reg.intercept_

输出：

(array([[0., 2., 5., 8.]]), array([-1.]))

算法设计2

# 固定numpy随机性
np.random.seed(42)
# 加上正态分布噪点
y = fx(x)+np.random.normal(-1,1,100).reshape(-1,1)
lin_reg = LinearRegression()#构建线性回归模型对象
lin_reg.fit(x_0,y)
lin_reg.coef_,lin_reg.intercept_

输出：

(array([[0.        , 2.14160723, 5.00509044, 7.99409475]]),
 array([-2.16616584]))

习题5.6

  函数$g(x)=\frac{10a}{10b+(a-10b)e^{-asinx}}$，取$a=1.1,b=0.01$，计算$x=1,2,\cdots ,20$时，$g(x)$对应的函数值，把这样得到的数据作为模拟观测值，记作$(x,y),i=1,2,\cdots ,20$，利用$i=1,2,\cdots ,20$反求解析式

算法设计

from scipy.optimize import curve_fit
def gx(x,a,b):
    return (10 * a)/(10 * b + (a - 10 * b)*np.exp(-a*np.sin(x)))
x = np.arange(1,21,1)
y = gx(x,1.1,0.01)
popt, pcov = curve_fit(gx, x, y)
popt

输出：

array([1.1 , 0.01])