如何用Sympy平滑绘制二阶ODE求解后的结构动力学响应曲线？

问题分析 #

你用Sympy求解的动力学方程响应曲线不够光滑，本质是Sympy默认的绘图采样点过少；而用Numpy时遇到的sin(5.3x)处理问题，可通过将Sympy符号表达式转换为Numpy可调用函数解决。

方案一：直接优化Sympy的plot函数 #

Sympy的plot()函数支持通过num_points参数增加采样点，或者开启adaptive自适应采样，来提升曲线光滑度：

from sympy import *
import matplotlib.pyplot as plt

x, y = symbols("x, y")
f = symbols('f', cls=Function)
y = f(x)

# 定义控制方程
eq = Eq(6938.78*diff(y, x, 2) + 5129.907*diff(y, x) + 379259*y - 900*sin(5.3*x), 0)
# 求解方程（带初始条件）
eq_done = dsolve(eq, y, ics={f(0): 0, diff(y, x).subs(x, 0): 0})

# 增加采样点至2000，开启自适应采样
plot(eq_done.rhs, (x, 0, 10), num_points=2000, adaptive=True)
plt.show()

方案二：转换为Numpy函数后用Matplotlib绘制（更灵活） #

用Sympy的lambdify工具将符号表达式转换为Numpy兼容的函数，再生成密集的x数组计算响应值，即可得到光滑曲线：

from sympy import *
import numpy as np
import matplotlib.pyplot as plt

x, y = symbols("x, y")
f = symbols('f', cls=Function)
y = f(x)

# 定义并求解方程
eq = Eq(6938.78*diff(y, x, 2) + 5129.907*diff(y, x) + 379259*y - 900*sin(5.3*x), 0)
eq_done = dsolve(eq, y, ics={f(0): 0, diff(y, x).subs(x, 0): 0})

# 将Sympy表达式转换为Numpy可调用函数
y_func = lambdify(x, eq_done.rhs, 'numpy')

# 生成密集的x数组（0到10，共1000个点）
x_vals = np.linspace(0, 10, 1000)
# 计算对应的y值
y_vals = y_func(x_vals)

# 绘制曲线
plt.figure(figsize=(10, 6))
plt.plot(x_vals, y_vals)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title('单自由度动力学方程响应曲线')
plt.grid(True)
plt.show()

这个方法完全解决了你提到的Numpy处理sin(5.3x)的问题，因为lambdify会自动将Sympy的sin转换为Numpy的np.sin，适配数值计算。

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