如何通过代码从已求解的ODE计算二阶ODE并修复位移求解异常

一、位移计算全零的问题及修复 #

问题原因 #

你用odeint计算位移时，velocity函数逻辑错误：odeint要求导数函数的第一个参数是当前待求解的状态量（此处为位移x），但你直接调用了外部变量u，实际运行时odeint传入的是当前x的数值（初始为0），导致dxdt始终为0，最终结果全零。

修复方案（三种可选） #

方案1：直接积分速度数组（最简单）

既然已得到速度u随时间t的变化数组，直接用累积积分计算位移：

from scipy.integrate import cumtrapz

# cumtrapz返回结果长度比t少1，补初始位移x0=0
x = np.concatenate([[x0], cumtrapz(u.flatten(), t)])

plt.xlabel('Time in s')
plt.ylabel('Displacement in m')
plt.plot(t, x)
plt.grid()
plt.show()

方案2：修正ODE函数，传入速度数组

若坚持用odeint，可通过插值获取当前时间对应的速度值：

def velocity(x, t, u, t_array):
    # 插值得到当前时间t的速度
    current_u = np.interp(t, t_array, u.flatten())
    return current_u

x = odeint(velocity, x0, t, args=(u, t))

方案3：联立求解速度与位移的ODE系统

将位移和速度作为状态向量，一次性求解更规范：

def system(state, t):
    x, u = state
    a = 3 * t**2  # 加速度表达式
    du_dt = a
    dx_dt = u
    return [dx_dt, du_dt]

# 初始状态：位移0，速度0
initial_state = [x0, u0]
result = odeint(system, initial_state, t)
x = result[:, 0]
u = result[:, 1]

plt.xlabel('Time in s')
plt.ylabel('Displacement in m')
plt.plot(t, x)
plt.grid()
plt.show()

二、Sympy导入错误的修复 #

你写的from sympy import symbol是拼写错误，正确导入为复数形式symbols：

# 导入单个符号
from sympy import symbols
# 或导入所有常用功能
from sympy import *

用Sympy做符号积分验证的示例代码：

import sympy as sp

t_sym = sp.symbols('t')
a_sym = 3 * t_sym**2  # 加速度符号表达式
u_sym = sp.integrate(a_sym, t_sym)  # 积分得速度
x_sym = sp.integrate(u_sym, t_sym)  # 积分得位移

# 转为数值函数代入t数组计算
x_func = sp.lambdify(t_sym, x_sym, 'numpy')
x = x_func(t)

plt.xlabel('Time in s')
plt.ylabel('Displacement in m')
plt.plot(t, x)
plt.grid()
plt.show()

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