《数学建模简明教程--基于python》学习笔记-第三章-插值与拟合-课后习题解答

本文链接：https://blog.csdn.net/delsystem32/article/details/123942435

文章目录

准备工作
1 分析插值与拟合的联系与区别
2 给出函数的n个节点，然后插值
3 已知四个数据点,然后插值求f(4)
4 非线性拟合求初始电压v0和充电常数 $\tau$
5 三次多项式拟合
6 非线性最小二乘法拟合
学习困惑

先纵览一下插值与拟合的题目
在这里插入图片描述

在这里插入图片描述

准备工作

导入绘图的库以及插值和拟合的函数。

import numpy as np
import matplotlib.pyplot as plt
# 插值
from scipy.interpolate import lagrange  #  拉格朗日插值
from scipy.interpolate import interp1d  #  一维插值:对一个自变量的函数进行插值

# 线性拟合 可以使用np.polyfit进行多项式拟合 (当然线性拟合不是只能是多项式只要参数线性，就可以用线性最小二乘拟合)
# 非线性拟合 使用curve_fit
from scipy.optimize import curve_fit
# 初始化设置
plt.rc('font', family='SimHei')     # 显示中文字体
plt.rc('axes',unicode_minus=False)  # 显示符号
np.set_printoptions(precision=3, suppress=True)   # 设置小数位数为3，并关闭print时的科学计数法
np.random.seed(2022)  # 设置随机数种子

1 分析插值与拟合的联系与区别

题1：
在这里插入图片描述

插值：给定一系列不连续的样本点，然后你需要根据这些点给出一种方法，它可以返回未知的更高密度的点，并且保留原来所有的点的位置不变（用你的插值方法）
拟合：用一种方法，它可以以最大解释程度来解释你所给出的点，可以是根据各种机理模型给出，也要有判断拟合好坏的标准，一般就是距离。
联系：都是要根据已知数据给出近似函数（插值函数和拟合函数）
区别：插值需要经过所有数据点，如用于地貌绘测等，本身在每个插值节点的数值往往是实实在在确定的、存在的，才会想到要一个经过所有数据点的函数，当然也会有观测误差存在。拟合则是按照一定的标准，再用某种方法来求得拟合效果最好的函数，所谓“拟合最好的函数”取决于拟合的标准，拟合较插值更有利于反映对象整体的变化趋势。

2 给出函数的n个节点，然后插值

使用拉格朗日，分段线性和三次样条三种插值方法

题2：
在这里插入图片描述

2.1 $y=sin(x),0\leqslant x\leqslant 2\pi$

# 已知数据点
x = np.append(np.random.random(7)*(np.pi*2),[0,np.pi*2])  # 我采用简单随机采样 但是添加了首位两个点 可以避免后续的一些插值范围问题
# 9个数据点都在[0,np.pi*2]上
# x = np.arange(0,np.pi*2+0.0001,np.pi/4)  # 等距采样 9个数据点
y = np.sin(x)

# 被插值点
x_ = np.arange(0,np.pi*2+0.0001,np.pi/20)

plt.figure(figsize=(14,4))

# 绘制拉格朗日8次插值，线性插值，三次样条插值，原始曲线sin(x)，以及采样点
y1 = lagrange(x, y)  # 插值8次多项式 返回的是每一个项数前的系数（从高到低）
plt.plot(x_, np.polyval(y1, x_), '--*b', label=f'拉格朗日{len(x)-1}次插值') 

f2 = interp1d(x, y); y2 = f2(x_) 
plt.plot(x_, y2, '-.', label='线性插值')  

f3 = interp1d(x, y, 'cubic'); y3 = f3(x_)
plt.plot(x_, y3, '-p', label='三次样条插值')

plt.plot(x_, np.sin(x_), 'r', label='原始曲线')

plt.scatter(x,y,c='k',s=100, label='数据点')

# 装饰性函数
plt.axis('equal')
plt.grid(axis='y')
plt.yticks(np.linspace(-1,1,5))
plt.xticks(np.linspace(0,np.pi*2,9))
plt.title('三种插值方法比较图')
plt.legend()
plt.show()

# 观察到：当取样点距离一样时，三次和拉格会重叠
# 在这个函数上 拉格朗日还算稳定的

在这里插入图片描述

增加插值节点
可以发现拉格朗日和三次样条的拟合曲线越来越接近
三种方法的拟合效果也越来越好
在插值节点较少时，拉格朗日插值会出现龙格现象
在插值节点数为9时，发现拉格朗日插值效果比三次样条插值更好，在一些有的两个点之间，三次样条的拟合曲线会更凸出去,但拉格朗日拟合的曲线却能较好逼近sin函数

2.2 $y=(1-x^2)^{1/2},-1\leqslant x\leqslant1$

# 已知点
x = np.append(np.random.random(5)*2-1,[-1,1])  # 在[0,np.pi*2]上随机采样，但是包括首位两个点，可以避免后续的一些插值范围问题
y = np.sqrt((1-x**2)) # 函数2

# 被插值点
x_ = np.linspace(-1,1)

plt.figure(figsize=(10,6))

# 绘制拉格朗日10次插值，线性插值，三次样条插值，原始曲线，以及采样点
y1 = lagrange(x, y)  # 插值10次多项式 返回的是每一个项数前的系数
plt.plot(x_, np.polyval(y1, x_), '--*b', label=f'拉格朗日{len(x)-1}次插值') 

f2 = interp1d(x, y); y2 = f2(x_) 
plt.plot(x_, y2, '-.', label='线性插值')

f3 = interp1d(x, y, 'cubic'); y3 = f3(x_)
plt.plot(x_, y3, '-p', label='三次样条插值')

plt.plot(x_, np.sqrt((1-x_**2)), 'r', label='原始曲线')

plt.scatter(x,y,c='k',s=100, label='数据点')

# 装饰性函数
plt.axis('equal')
plt.grid(axis='y')
plt.yticks(np.linspace(-0.2,1,10))
plt.xticks(np.linspace(min(x_),max(x_),10))
plt.title('三种插值方法比较图')
plt.legend()
plt.show()

#  观察到：拉格朗日有时会突出去 三次样条的效果还是不错的

在这里插入图片描述

2.3 $\cos^{10}x,-4\leqslant x\leqslant4$

def fun3(x):
    return np.cos(x)**10

# 这里我选择[-4,4]为定义域，主要想看看对周期函数插值的效果
x = np.append(np.random.random(8)*8-4,[-4,4])  # 在[0,np.pi*2]上随机采样，但是包括首位两个点，可以避免后续的一些插值范围问题
y = fun3(x)

# 被插值点
x_ = np.linspace(-4,4,100)
y_ = fun3(x_)

plt.figure(figsize=(10,6))

# 绘制拉格朗日9次插值，线性插值，三次样条插值，原始曲线，以及采样点
y1 = lagrange(x, y)  # 插值9次多项式 返回的是每一个项数前的系数
plt.plot(x_, np.polyval(y1, x_), '--*b', label=f'拉格朗日{len(x)-1}次插值') 

f2 = interp1d(x, y); y2 = f2(x_) 
plt.plot(x_, y2, '-.', label='线性插值')

f3 = interp1d(x, y, 'cubic'); y3 = f3(x_)
plt.plot(x_, y3, '-p', label='三次样条插值')

plt.plot(x_, y_, 'r', label='原始曲线')

plt.scatter(x,y,c='k',s=100, label='已知数据点')

# 装饰性函数
plt.axis('equal')
plt.grid(axis='y')
plt.yticks(np.linspace(-0.2,1,10))
plt.xticks(np.linspace(min(x_),max(x_),10))
plt.ylim(min(y_)*1.5, max(y_)*1.5)
plt.title('三种插值方法比较图')
plt.legend()
plt.show()

#  观察到：拉格朗日有时会凸出去 三次样条的效果也勉强

在这里插入图片描述

2.4 $y=e^{-x^2},-2\leqslant x\leqslant2$

def fun4(x):
    return np.exp(-x**2)

# 已知点
x = np.append(np.random.random(5)*4-2,[-2,2])  # 在[0,np.pi*2]上随机采样，但是包括首位两个点，可以避免后续的一些插值范围问题
y = fun4(x)

# 被插值点
x_ = np.linspace(-2,2,50)
y_ = fun4(x_)

p = (max(y_)-min(y_))/(max(x_)-min(x_))
plt.figure(figsize=(10,p*10))

# 绘制拉格朗日6次插值，线性插值，三次样条插值，原始曲线，以及采样点
y1 = lagrange(x, y)  # 插值10次多项式 返回的是每一个项数前的系数
plt.plot(x_, np.polyval(y1, x_), '--*b', label=f'拉格朗日{len(x)-1}次插值') 

f2 = interp1d(x, y); y2 = f2(x_) 
plt.plot(x_, y2, '-.', label='线性插值')

f3 = interp1d(x, y, 'cubic'); y3 = f3(x_)
plt.plot(x_, y3, '-p', label='三次样条插值')

plt.plot(x_, y_, 'r', label='原始曲线')

plt.scatter(x,y,c='k',s=100, label='数据点')

# 装饰性函数
plt.axis('equal')
plt.grid(axis='y')
plt.yticks(np.linspace(-0.2,1,7))
plt.xticks(np.linspace(min(x_),max(x_),10))
plt.ylim(-0.2, max(y_)*1.1)
plt.title('三种插值方法比较图')
plt.legend()
plt.show()

# 观察到：拉格朗日有时会突出去 三次样条的效果还是不错的

在这里插入图片描述

3 已知四个数据点,然后插值求f(4)

题3：
在这里插入图片描述

# 生成已知数据
x = [1,3,5,7]
y = [1,8,10,3]
# 生成插值点
x_ = np.linspace(1,7)

a = lagrange(x, y)  # 拉格朗日插值三次多项式 a 是高次到低次的

f = interp1d(x, y) # 线性插值 返回的是一个插值函数

g = interp1d(x, y, 'cubic')  # 三次样条插值

# 绘图
plt.plot(x_, np.polyval(a, x_), '--*b', label='拉格朗日3次插值')
plt.plot(x_, f(x_), '-.', label='线性插值')
plt.scatter(x_, g(x_), c='None', edgecolors='r', s=100,label='三次样条插值')
plt.scatter(x, y, c='r', s=100, label='已知数据点')
plt.legend()
plt.show()

在这里插入图片描述

print(f'用三种方法试求f(4)\n    线性插值:{np.polyval(a, 4)}\n    拉格朗日3次插值:{f(4)}\n    三次样条插值:{g(4)}')

用三种方法试求f(4)
    线性插值:9.874999999999996
    拉格朗日3次插值:9.0
    三次样条插值:9.875000000000002

4 非线性拟合求初始电压v0和充电常数 $\tau$

题4：
在这里插入图片描述

x0 = np.array([0.5,1,2,3,4,5,7,9])
y0 = np.array([6.36,6.48,7.26,8.22,8.66,8.99,9.43,9.63])
p = np.polyfit(x0, y0, 2)

def fun0(t,v0,r):
        return 10 - (10 - v0)*np.exp(-t/r)

from scipy.optimize import curve_fit  # curve_fit的详细内容在题6介绍

popt, pcov = curve_fit(fun0, x0, y0)
popt

array([5.558, 3.5  ])

解得v0=5.558, $\tau$ =3.5

5 三次多项式拟合

题5：
在这里插入图片描述

#  生成数据
x = np.linspace(-1, 1, 11)
y = np.array([1, 1.08, 0.98, 0.73, 0.39, 0 ,-0.39, -0.73, -0.98, -1.08, -1])

plt.plot(x, y, 'o')
plt.show()

先看一下数据的分布情况。
在这里插入图片描述

这里我经过多次试验，然后决定用三次多项式来拟合

# 生成100个节点用于绘制拟合曲线
x_ = np.linspace(-3, 3, 100) 
# 三次多项式拟合 
p = np.polyfit(x,y,3)  
Y = np.polyval(p, [-0.85, 0.1, 0.85])  # 求题目要求的目标值
# 绘制拟合图像，展示拟合效果
plt.plot(x_, np.polyval(p, x_), label = '拟合曲线')
plt.plot(x, y, 'o', label = '已知数据')
plt.scatter([-0.85, 0.1, 0.85], Y, c='g', s=100, label='目标值')
plt.title('已知数据点和三次拟合曲线')
plt.axis('equal')
plt.ylim(-2,2)
plt.legend()
# 添加数据标签
plt.text(-0.85, Y[0]+0.2, f'(-0.85,{np.round(Y[0],2)})', c='g')
plt.text(0, Y[1]+0.2, f'(-0.1,{np.round(Y[1],2)})', c='g')
plt.text(0.85, Y[2]+0.2, f'(0.85,{np.round(Y[2],2)})', c='g')
plt.show()

在这里插入图片描述

print('x=-0.85,0.1,0.85时的函数值:',np.polyval(p, [-0.85, 0.1, 0.85]))

x=-0.85,0.1,0.85时的函数值: [ 1.082 -0.197 -1.082]

6 非线性最小二乘法拟合

题6：
在这里插入图片描述

6.1 生成数据并定义拟合函数/模型

# 生成已知数据
x = np.linspace(0.1,1,10)
y = np.array([2.32, 2.64, 2.97, 3.28, 3.6, 3.9, 4.21, 4.51, 4.82, 5.12])

# 定义已知的拟合函数/模型
def f(x, a, b, c, d):  # 注意函数第一个参数在拟合时会被当作自变量，后面的才是要拟合的参数
    return a*x + b*(x**2)*np.exp(-c*x) + d  # 参数非线性

6.2 非线性最小二乘法求解参数

popt, pcov = curve_fit(f, x, y)

popt 是拟合的参数 pcov 是参数的协方差矩阵
官方文档有写Default is ‘lm’ for unconstrained problems and ‘trf’ if bounds are provided.
也就是说在参数没有限制定义域时，会默认使用lm方法，且参数初值默认为1，在有参数范围限制时为trf方法

点击查看 curve_fit 官方文档

print('a, b, c, d的拟合值为:', np.round(popt, 4))

a, b, c, d的拟合值为: [ 3.397 -0.464  0.582  1.982]

6.3 绘制拟合曲线

x_ = np.linspace(min(x), max(x))  # 生成数据点来画拟合函数
y_ = f(x_, *popt)  # *popt表示参数传递 会把popt中的值作为参数一个个传入 

plt.plot(x, y, 'o', label= '已知数据')
plt.plot(x_, y_, label='拟合曲线')

plt.title('已知数据的散点图和拟合的曲线')
plt.axis('square')
plt.legend()
plt.show()