基物实验 Python 曲线拟合的应用

问题背景

作者在做基物实验 波尔共振仪研究受迫振动 的实验中，由于其操作不当，或者读数时没长眼睛，导致得到的幅频特性曲线呈现如下惨状：

这么惨不忍睹的曲线，不加处理地交上去肯定是得不到合格的分数的……

于是他决定进行曲线拟合。

---

解决方法

1. 确定待拟合函数模式

首先，我所掌握的方法中，一般都需要我们清楚要拟合的曲线的基本表达式。

根据预习报告中的相关知识，我们得到有用的信息：
$$\theta =\frac{h}{\sqrt{({\omega }_0^2-{\omega }^2{)}^2+4{\beta }^2{\omega }^2}}$$
在幅频特性曲线中，横坐标 $x$ 代表 $\omega /{\omega }_0$ ，因此表达式可以转化成：
$$\theta =\frac{h}{\sqrt{{\omega }_0^4(1-x^2{)}^2+4{\beta }^2{\omega }_0^2{x}^2}}$$
${\omega }_0$ 是实验的前几个步骤中所测定的常数，所以待定参数只有两个——$h$ 和 $b:=4{\beta }^2$ 。

所以我们得到了函数模板：

from numpy import *

def function(x, b, h):
    result = [h/(sqrt(w0**4 * (1-xi**2)**2 + b * w0**2 * xi**2)) for xi in x]
    return array(result)

2. 最小二乘拟合

使用 scipy.optimize 中的 leastsq 函数进行曲线最小二乘拟合。

$leastsq$ 实际上是一个最小器，用来求函数的最小值的。而被求最小值的函数即为我们的误差函数。

def ferr(p, func, y, x): #误差函数，p参数包
    return y - func(x, *p)

下面我们要做的就是根据已有数据列表求出使得误差函数最小的待定参数。

def myfit(func, n, x, y):  #func要拟合的函数模板，n待定参数个数
    x, y = array(x), array(y)
    c, ret_val = lsq(ferr, [0]*n, args=(func, y, x))
    print(ret_val)
    if ret_val in range(1, 5):
        return c  #返回合适的参数包

$ret\_val$ 取值在 $0\sim 4$ 之间为拟合成功。

这里我把所有参数的初始值设成了 $0$ ，具体迭代初始值设为多少视情况而定。

返回的 $c$ 即合适的参数包， 只需在形参中展开就能得到我们想要的函数了。

3. 拟合验证

根据物理学相关知识，$\omega /{\omega }_0$ 取 $1$ 时，我们称这种状态为 共振 ，此时角幅值 $\theta$ 达到最大值。

那么下面我们只需要验证拟合得得函数的最大值是否在 $1$ 处取得即可。

scipy.optimize 中提供了一款非常简便得求某点附近极值的函数 fmin。

但是很不幸，scipy 没有提供求最大值的函数，所以这时就需要我们手动改写一下 fmin 了。

参数照 fmin 函数原方不动照搬即可，或者可以加入我们给函数赋予待定参数包的形参接口：

def fmax(func, x0, args=(), **kwargs):
    return fmin(lambda x: -func(x, *args), x0, **kwargs)  #加个负号相当于求极大值

求得取最大值的位置分别是 $0.99958496$ 和 $0.99956055$ ，和 $1$ 非常接近，可见拟合精度非常之高。

---

结果

完成了，可见已知函数形式的曲线拟合非常简单易行。

在自己实验数据测得不如人家的时候，可以使用如下办法扳回一局。

附上全部代码以及拟合效果：

from  numpy import *
import math
import matplotlib.pyplot as plt
from scipy.optimize import leastsq as lsq

plt.rcParams['font.sans-serif']=['SimHei']  #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False  #用来正常显示负号

# 最小二乘拟合
w0 = 2*pi/1.5499
def function(x, b, h):
    result = [h/(sqrt(w0**4 * (1-xi**2)**2 + b * w0**2 * xi**2)) for xi in x]
    return array(result)

def ferr(p, func, y, x): #误差函数，p参数包
    return y - func(x, *p)

def myfit(func, n, x, y):  #func要拟合的函数模板，n待定参数个数
    x, y = array(x), array(y)
    c, ret_val = lsq(ferr, [0]*n, args=(func, y, x))
    print(ret_val)
    if ret_val in range(1, 5):
        return c  #返回合适的参数包
# 读入数据
fp = open("data.in")
phi1 = list(map(float, fp.readline().split()))
T1 = list(map(float, fp.readline().split()))
phi2 = list(map(float, fp.readline().split()))
T2 = list(map(float, fp.readline().split()))
n = len(T1) 
w1 = [2*pi/Ti/w0 for Ti in T1]
w2 = [2*pi/Ti/w0 for Ti in T2]

# 原图展示
# plt.subplot(1, 2, 1)
# plt.plot(w1, phi1, linestyle='--', marker='o')
# plt.plot(w2, phi2, linestyle='-.', marker='^')

# plt.subplot(1, 2, 2)
plt.plot(w1, phi1, 'b*', label = r"阻尼1")
new_w1 = linspace(start = 0.95, stop = 1.05, num = 25)
para1 = myfit(function, 2, w1, phi1)
plt.plot(new_w1, function(new_w1, *para1), 'g--', label = r"阻尼1拟合曲线")

plt.plot(w2, phi2, 'rx', label = r"阻尼3")
new_w2 = linspace(start = 0.95, stop = 1.05, num = 25)
para2 = myfit(function, 2, w2, phi2)
plt.plot(new_w2, function(new_w2, *para2), 'k-', label = r"阻尼3拟合曲线")

# 求最大值
from scipy.optimize import fmin

def fmax(func, x0, args=(), **kwargs):
    return fmin(lambda x: -func(x, *args), x0, **kwargs)  #加个负号相当于求极大值

result = fmax(function, 1, para1)
print(result)

result = fmax(function, 1, para2)
print(result)

plt.grid()
plt.legend()
plt.show()

效果对比展示：

验证结果：