多次使用scipy.solve_ivp求解变扭矩ODE异常问题排查

TLDR：使用scipy.solve_ivp求解电机驱动旋转机翼扭矩的ODE，结果数值异常（phi'达1e12、phi''达1e26），无法定位原因。

我需要求解的方程为：phi'' = a*torque_{motor}-b*(phi')^2，其中a、b是常数，torque_{motor}会在每次迭代中更新。我的思路是分时间窗口[t,t+dt]、[t+dt,t+2dt]...循环求解该ODE，每个窗口内给torque_motor赋予一个新的浮点值，以此模拟电机在不同时间步施加不同力矩。

最小可运行代码 #

import numpy as np
from scipy.integrate import solve_ivp

class Solver():
    def __init__(self):
        self.torque = 0

    def set_torque(self, new):
        self.torque = new

    def phi_derivatives(self, t, y):
        """
        定义待求解的ODE：I * phi_ddot = tau_z - tau_drag.
        将y视为向量 y = [phi,phi_dot]，求解 dy/dt = f(y,t)
        """
        # 常数取值
        a = 1090179.6616082331
        b = 271.9406850484702

        phi, phi_dot = y[0], y[1]
        dy_dt = [phi_dot, a * self.torque - b * (phi_dot ** 2)]
        return dy_dt

# 初始化存储数组
phi_arr = []
phi_dot_arr = []
phi_ddot_arr = []
phi_0 = 0
phi_dot_0 = 0.01
start_t = 0
end_t = 0.05
delta_t = 0.001
s = Solver()
# 生成正弦波扭矩序列
sin_torques = np.sin(2 * np.pi * np.linspace(0, 1, 20))

for torque in sin_torques:
    s.set_torque(torque)
    sol = solve_ivp(s.phi_derivatives, t_span=(start_t, end_t), y0=[phi_0, phi_dot_0])
    phi, phi_dot = sol.y
    _, phi_ddot = s.phi_derivatives(0, [phi, phi_dot])
    phi_arr.append(phi)
    phi_dot_arr.append(phi_dot)
    phi_ddot_arr.append(phi_ddot)

    # 更新下一个时间窗口的初始条件和时间范围
    phi_0 = phi[-1]
    phi_dot_0 = phi_dot[-1]
    start_t = end_t
    end_t += 0.05

# 合并数组
phi_arr = np.concatenate(phi_arr)
phi_dot_arr = np.concatenate(phi_dot_arr)
phi_ddot_arr = np.concatenate(phi_ddot_arr)

问题 #

每个时间窗口的扭矩值取自正弦波，我原本预期phi会随扭矩方向变化产生振荡，但实际得到的phi'和phi''数值异常，量级分别达到1e12和1e26。请问这是实现上的问题吗？

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