如何用SymPy绘制dsolve求解的微分方程解的函数图像

以下是对应Matlab代码的Python实现，基于SymPy完成符号计算、绘图及求解：

步骤1：导入所需库 #

import sympy as sp
import matplotlib.pyplot as plt
import numpy as np

步骤2：定义符号变量与函数 #

对应Matlab的syms y(t) r V：

t = sp.symbols('t')
r, V = sp.symbols('r V')
y = sp.Function('y')(t)

步骤3：构建微分方程与初始条件 #

对应Matlab的eqn = diff(y) == -r / V * y和cond = y(0) == 1：

eqn = sp.Eq(sp.diff(y, t), -r / V * y)
cond = sp.Eq(y.subs(t, 0), 1)

步骤4：求解微分方程 #

对应Matlab的y = dsolve(eqn, cond)：

sol = sp.dsolve(eqn, cond)
y_sol = sol.rhs  # 提取解的右侧表达式

步骤5：代入参数数值 #

对应Matlab的两次subs操作：

r_val = 3.663959132E10
V_val = 4871E9
y_sub = y_sol.subs([(r, r_val), (V, V_val)])

步骤6：绘制函数图像 #

对应Matlab的fplot(y, [1, 1000])，这里采用更灵活的Matplotlib绘制方式：

# 将符号函数转换为可计算的数值函数
y_numeric = sp.lambdify(t, y_sub, 'numpy')
t_vals = np.linspace(1, 1000, 1000)
y_vals = y_numeric(t_vals)

plt.plot(t_vals, y_vals)
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('Solution of Differential Equation')
plt.show()

步骤7：求解y=0.05时的t值 #

对应Matlab的val = double(solve(y == 0.05, t))：

t_sol = sp.solve(sp.Eq(y_sub, 0.05), t)[0]
val = float(t_sol.evalf())  # 转换为高精度浮点数
print(f"y=0.05时的t值：{val}")

完整整合代码 #

import sympy as sp
import matplotlib.pyplot as plt
import numpy as np

# 定义符号变量与函数
t = sp.symbols('t')
r, V = sp.symbols('r V')
y = sp.Function('y')(t)

# 构建微分方程与初始条件
eqn = sp.Eq(sp.diff(y, t), -r / V * y)
cond = sp.Eq(y.subs(t, 0), 1)

# 求解微分方程
sol = sp.dsolve(eqn, cond)
y_sol = sol.rhs

# 代入参数
r_val = 3.663959132E10
V_val = 4871E9
y_sub = y_sol.subs([(r, r_val), (V, V_val)])

# 绘图
y_numeric = sp.lambdify(t, y_sub, 'numpy')
t_vals = np.linspace(1, 1000, 1000)
y_vals = y_numeric(t_vals)

plt.plot(t_vals, y_vals)
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('Solution of Differential Equation')
plt.show()

# 求解y=0.05时的t值
t_sol = sp.solve(sp.Eq(y_sub, 0.05), t)[0]
val = float(t_sol.evalf())
print(f"y=0.05时的t值：{val}")

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