Sympy lambdify转换符号表达式报错，求解ODE绘制吸引域

错误根源 #

1. 未定义符号变量：错误表达式中存在未处理的符号s（以及可能的d、ep），这些符号未被赋值或替换，导致lambdify无法生成合法的数值函数。
2. 错误的函数替换方式：使用math.sin/math.cos替换SymPy的三角函数，而math模块函数仅支持单个数值输入，无法处理solve_ivp传入的NumPy数组，引发类型错误。
3. ODE函数签名不匹配：solve_ivp要求导数函数的签名为func(t, y)，若你的ODE是dy/dx = v1(x,y)，则需要调整函数参数对应关系。

分步解决代码 #

1. 清理符号表达式，替换未知常数 #

首先确保所有符号变量（如s、d、ep）都被替换为具体数值：

import sympy as sp
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义所有符号变量
x, y, s, d, ep = sp.symbols('x y s d ep')

# 替换所有未知常数为实际数值（根据你的模型修改）
v1 = v1.subs(s, 1.0)
v1 = v1.subs(d, 0.5)
v1 = v1.subs(ep, 0.1)

# 简化表达式（可选，提升计算效率）
v1 = sp.simplify(v1)

2. 正确使用lambdify生成数值函数 #

指定modules='numpy'，让SymPy自动将三角函数转换为支持数组运算的NumPy版本：

# 生成适配NumPy的数值函数
v1_func = sp.lambdify((x, y), v1, modules='numpy')

# 调整函数签名适配solve_ivp：输入为(t, y)，输出为dy/dt
def ode_func(t, y):
    # t对应x，y[0]是当前y值
    return v1_func(t, y[0])

3. 求解ODE并绘制结果 #

# 求解参数设置
x_vals = np.linspace(0, 10, 100)
initial_y = [0.0]

# 求解ODE
solution = solve_ivp(ode_func, (x_vals[0], x_vals[-1]), initial_y, t_eval=x_vals)

# 绘制轨迹
plt.plot(x_vals, solution.y[0], 'r-', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.title('ODE Solution Trajectory')
plt.grid(True)
plt.show()

4. 绘制Lyapunov函数的吸引域 #

通过遍历初始条件，绘制向量场与轨迹：

# 定义初始条件范围
x_min, x_max = -5, 5
y_min, y_max = -5, 5
x_grid = np.linspace(x_min, x_max, 25)
y_grid = np.linspace(y_min, y_max, 25)
X, Y = np.meshgrid(x_grid, y_grid)

# 计算向量场（dx/dt=1，dy/dt=v1(x,y)）
U = np.ones_like(X)
V = np.zeros_like(Y)
for i in range(X.shape[0]):
    for j in range(X.shape[1]):
        V[i,j] = v1_func(X[i,j], Y[i,j])

# 绘制向量场
plt.quiver(X, Y, U, V, color='lightgray', alpha=0.6)

# 绘制多个初始条件的轨迹
for y0 in np.linspace(y_min, y_max, 8):
    sol = solve_ivp(ode_func, (x_min, x_max), [y0], t_eval=np.linspace(x_min, x_max, 50))
    plt.plot(sol.t, sol.y[0], 'b-', alpha=0.7)

# 标记平衡点（替换为你模型的实际平衡点）
plt.scatter(0, 0, color='red', marker='*', s=150, label='Equilibrium')

plt.xlabel('x')
plt.ylabel('y')
plt.title('Lyapunov Function Attraction Domain')
plt.legend()
plt.grid(True)
plt.show()

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