使用scipy的solve_bvp求解含加速度边值问题遇奇异雅可比矩阵错误

问题 #

用scipy.integrate.solve_bvp从加速度数据推导速度和位移时，设置初末时刻速度均为0的边界条件会触发A singular Jacobian encountered when solving the collocation system错误；但将边界条件改为初始和末时刻位移均为0时，代码可正常运行。

报错原因 #

咱们的微分方程组是：

位移的导数 = 速度
速度的导数 = 加速度

如果只约束初末速度为0，位移就存在“自由度”——所有位移值加任意常数，速度依然满足条件（因为导数不变）。这种解不唯一的情况会导致求解器的雅可比矩阵奇异，无法计算出唯一解。

而约束初末位移为0时，位移的解被固定，速度也随之确定，求解器就能正常收敛。

解决办法 #

如果实际场景需要初末速度为0，得给位移加一个约束（比如固定初始位移为0）来确保解的唯一性；同时不要用全零矩阵当初始猜测，改用累积积分的结果作为初始值，帮助求解器更快收敛：

1. 修改边界条件：固定初始位移s(0)=0 + 末时刻速度v(T)=0；
2. 用累积积分得到的位移、速度替换全零矩阵作为初始猜测。

修改后的代码 #

from scipy.integrate import solve_bvp, cumulative_trapezoid
from scipy.interpolate import UnivariateSpline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# 读取加速度数据并预处理
data = pd.read_csv("Linear Accelerometer.csv", header=None, skiprows=1)
t_raw = data.iloc[:, 0].astype(str).str.replace(',', '.').astype(float).values
a_raw = data.iloc[:, 2].astype(str).str.replace(',', '.').astype(float).values

n_points = 1000
t_coarse = np.linspace(t_raw[0], t_raw[-1], n_points)
a_coarse = np.interp(t_coarse, t_raw, a_raw)

# 平滑加速度曲线
a_spline = UnivariateSpline(t_coarse, a_coarse, s=0.1)
a_func = a_spline

# 生成初始猜测：用累积积分得到位移和速度
v_guess = cumulative_trapezoid(a_func(t_coarse), t_coarse, initial=0)
s_guess = cumulative_trapezoid(v_guess, t_coarse, initial=0)
y_init = np.vstack([s_guess, v_guess])  # 替换全零矩阵

# 定义微分方程组
def fun(t, y):
    return np.array([y[1], a_func(t)])

# 新边界条件：初始位移为0，末时刻速度为0
def bc(ya, yb):
    return np.array([ya[0], yb[1]])

# 调用求解器
sol = solve_bvp(fun, bc, t_coarse, y_init, max_nodes=100000, verbose=2)

# 打印求解状态
print(sol.message)

# 提取结果
s = sol.sol(t_coarse)[0]
v = sol.sol(t_coarse)[1]

# 绘图展示
plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(t_coarse, s, label="位移 s(t)")
plt.ylabel("位移 [m]")
plt.legend()
plt.grid()

plt.subplot(2, 1, 2)
plt.plot(t_coarse, v, label="速度 v(t)", color="orange")
plt.xlabel("时间 [s]")
plt.ylabel("速度 [m/s]")
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

补充说明 #

如果必须严格约束初末速度都为0，还可以添加位移的相对约束（比如s(T) - s(0) = 0），但需要结合实际物理场景判断合理性。通常固定初始位移为0是更贴合实际的选择。

内容的提问来源于stack exchange，提问作者Luis Utermann