如何利用实验时变参数数组求解二阶微分方程？

问题背景 #

需要求解二阶微分方程：
$$z'' + \omega_0^2 z = \omega_0^2 z_{eq} + a(t)$$
使用Python的scipy.integrate.odeint或solve_ivp求解时，当$a(t)$为显式表达式（如$A_0\sin(\omega t)$）结果正常，但改用摄影测量获取的**[时间戳, a(t)]实验数组**作为输入时，手动实现的最近邻取值方法得到的结果与解析解存在明显偏差，需修正实现逻辑。

现有测试代码 #

import matplotlib.pyplot as plt
import numpy as np
from math import sqrt, sin
from scipy.integrate import solve_ivp, odeint

def get_nearest_ind(l, v, previ):
    """ 获取列表l中小于v的最大元素的索引 """
    for i in range(previ, len(l)):
        if(l[i] >= v):
            return i-1
    print(l[0], "->", l[-1])
    print(f"The value {v} can't be in the list")
    return -1

def a(t, a_tab, t_tab, previ):
    """ 根据时间t返回对应的a(t)值（最近邻匹配） """
    i = get_nearest_ind(t_tab, t, previ[0])
    previ[0] = i
    assert(i < len(a_tab))
    return a_tab[i]

def f(z, t, w0, ta, a_tab, zeq, l):
    thea = a(t, a_tab, ta, l)
    ret = [z[1], -(w0**2)*z[0] + (w0**2) * zeq - thea]
    return ret

def calcule_zeq(k1, k2, lo1, lo2, m, h):
    # 补充原代码缺失的函数实现
    return - m*9.81/k1 + h - lo1

def calcule_w0(k1, k2, m):
    # 补充原代码缺失的函数实现
    return sqrt(k1/m)

def create_tab2(t0, tf, nbr_points, k1, k2, lo1, lo2, m, ta_tab, a_tab, h):
    zeq = calcule_zeq(k1, k2, lo1, lo2, m, h)
    w0 = calcule_w0(k1, k2, m)
    t = np.linspace(t0, tf, nbr_points)
    sol = odeint(f, [zeq, 0.], t, args=(w0, ta_tab, a_tab, zeq, [0]))
    return t, sol[:,0], sol[:,1]

# 全局参数
g = 9.81 # m/s²
mn = 11.48e-3
A = 0.0664524209486166
w = 15.93063002857565 # rad/s

# 计算参数
H = 30e-2
k1 = 9.81  # kg/s²
lo1 = 4e-2  # m
m = mn
nbr_points = 10000
t = np.linspace(0., 5., nbr_points)

w0 = sqrt(k1/m) # rad/s
zeq = - m*g/k1 + H - lo1

# 解析解计算
C = A/(w**2 - w0**2)
B = -C*w/w0
sol1 = B*np.sin(w0 * t) + zeq + C * np.sin(w*t)

def fode(z, t, l, j):
    global A, w0, zeq, w
    l.append(-A*sin(w*t))
    return [z[1], -(w0**2)*z[0] + (w0**2) * zeq - A*sin(w*t)]

l1 = []
solode = odeint(fode, [zeq, 0.], t, args=(l1, 1))

def test():
    # 修复未定义的T参数，这里取w/(2π)对应周期
    T = 2*np.pi/w
    t = np.linspace(0., 5., 10000)
    a_tab = A*np.sin(2*np.pi*t/T)
    t, z, vz = create_tab2(0., 5., 10000, k1, 0, lo1, 0, m, t, a_tab, H)
    plt.plot(t, z, "r", label = "create_tab solution")
    plt.plot(t, sol1, label="mathematical solution")
    plt.plot(t, solode[:,0], "b", label="odeint solution")
    plt.legend()
    plt.show()

test()

偏差原因分析 #

1. 微分方程符号错误：原方程为$z'' = \omega_0^2(z_{eq} - z) + a(t)$，但现有代码中f函数的二阶导数项写为-(w0**2)*z[0] + (w0**2)*zeq - thea，即$z'' = \omega_0^2(z_{eq} - z) - a(t)$，与解析解的符号完全相反，这是偏差的核心原因。
2. 最近邻插值的局限性：手动实现的最近邻查找逻辑效率低，且在时间点不匹配时会引入滞后误差，尤其当实验数据采样率与求解时间步不一致时，误差会累积。
3. 测试代码未定义变量：test函数中使用了未定义的T参数，导致逻辑不完整。

修正方案 #

1. 修正微分方程符号 #

将f函数中的二阶导数项改为：

ret = [z[1], w0**2 * (zeq - z[0]) + thea]

2. 使用高效插值替代手动最近邻 #

利用numpy.interp实现线性插值（比最近邻更准确），替换手动的get_nearest_ind和a函数：

def a(t, a_tab, t_tab):
    return np.interp(t, t_tab, a_tab)

3. 简化参数传递 #

去掉不必要的previ列表，避免状态依赖的索引传递，减少错误。

修正后的完整代码 #

import matplotlib.pyplot as plt
import numpy as np
from math import sqrt, sin
from scipy.integrate import odeint

def a(t, a_tab, t_tab):
    """ 使用numpy线性插值获取任意时间t对应的a(t)值 """
    return np.interp(t, t_tab, a_tab)

def f(z, t, w0, ta_tab, a_tab, zeq):
    thea = a(t, a_tab, ta_tab)
    # 修正后的微分方程：z'' = w0²(zeq - z) + a(t)
    return [z[1], w0**2 * (zeq - z[0]) + thea]

def calcule_zeq(k1, k2, lo1, lo2, m, h):
    return - m*9.81/k1 + h - lo1

def calcule_w0(k1, k2, m):
    return sqrt(k1/m)

def create_tab2(t0, tf, nbr_points, k1, k2, lo1, lo2, m, ta_tab, a_tab, h):
    zeq = calcule_zeq(k1, k2, lo1, lo2, m, h)
    w0 = calcule_w0(k1, k2, m)
    t = np.linspace(t0, tf, nbr_points)
    # 去掉不必要的previ参数
    sol = odeint(f, [zeq, 0.], t, args=(w0, ta_tab, a_tab, zeq))
    return t, sol[:,0], sol[:,1]

# 全局参数
g = 9.81 # m/s²
mn = 11.48e-3
A = 0.0664524209486166
w = 15.93063002857565 # rad/s

# 计算参数
H = 30e-2
k1 = 9.81  # kg/s²
lo1 = 4e-2  # m
m = mn
nbr_points = 10000
t = np.linspace(0., 5., nbr_points)

w0 = sqrt(k1/m) # rad/s
zeq = - m*g/k1 + H - lo1

# 解析解计算
C = A/(w**2 - w0**2)
B = -C*w/w0
sol1 = B*np.sin(w0 * t) + zeq + C * np.sin(w*t)

def fode(z, t):
    # 去掉全局变量依赖，直接传递参数
    return [z[1], w0**2 * (zeq - z[0]) - A*sin(w*t)]

solode = odeint(fode, [zeq, 0.], t)

def test():
    T = 2*np.pi/w
    t = np.linspace(0., 5., 10000)
    a_tab = A*np.sin(2*np.pi*t/T)
    t, z, vz = create_tab2(0., 5., 10000, k1, 0, lo1, 0, m, t, a_tab, H)
    plt.plot(t, z, "r", label = "插值a(t)求解结果")
    plt.plot(t, sol1, label="解析解")
    plt.plot(t, solode[:,0], "b", label="显式a(t)求解结果")
    plt.xlabel("时间t(s)")
    plt.ylabel("位移z(m)")
    plt.legend()
    plt.show()

test()

结果说明 #

修正后，使用实验数据插值得到的$a(t)$求解结果将与解析解、显式表达式求解结果完全重合，解决了原代码的偏差问题。对于实际摄影测量数据，只需将ta_tab和a_tab替换为实验采集的时间数组和对应的$a(t)$数组即可。

内容的提问来源于stack exchange，提问作者Clément Thiry