如何用Matplotlib绘制Sympy隐函数并解决线条扩散问题

需求背景 #

我有一个隐函数（比如x**2 - y = 0），想在指定x范围内绘制图像，但sympy.plot_implicit生成的线条存在扩散问题，效果不满意。我希望获取绘图的数值数据，因此更倾向使用pyplot.plot。我能熟练用以下代码绘制显式SymPy函数，但不知道如何处理exp = sym.Eq(x**2 - y, 0)这类隐函数：

import sympy as sym
import numpy as np
from matplotlib import pyplot as plt

x, y = sym.symbols('x y', nonnegative=True)
exp = x**2

# 转换为NumPy可用函数绘图
x_arr = np.linspace(-2, 2, 100)
exp_func = sym.lambdify(x, exp, 'numpy')
exp_arr = exp_func(x_arr)

plt.plot(x_arr, exp_arr)

实际复杂表达式问题 #

我的实际需求是绘制方程b_sim = -1的图像，其中b_sim表达式如下：

from sympy import *

h, nu = symbols('h nu', nonnegative=True) 
b_sim = 1.0*cos(pi*sqrt(1 - h)/(2*nu))*cos(pi*sqrt(h + 1)/(2*nu)) - 1.0*sin(pi*sqrt(1 - h)/(2*nu))*sin(pi*sqrt(h + 1)/(2*nu))/sqrt(1 - h**2)

用sym.plot_implicit(b_sim + 1, (nu,0.225,1.5), (h, -1.1, 1.1))绘图时，线条扩散问题很明显。尝试用roots函数求解方程Eq(b_sim + 1, 0)的解析解，但直接报错。

一、简化隐函数的处理（可解析求解） #

对于像x² - y = 0这类能直接解出显式表达式的隐函数，步骤如下：

1. 用SymPy求解隐函数中目标变量关于另一变量的表达式
2. 将解析解转换为NumPy可用函数
3. 生成数值数组后用pyplot.plot绘图

示例代码：

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

x, y = sym.symbols('x y')
# 定义隐函数方程
eq = sym.Eq(x**2 - y, 0)
# 求解y关于x的表达式
sol_y = sym.solve(eq, y)[0]

# 转换为NumPy函数
x_arr = np.linspace(-2, 2, 100)
y_func = sym.lambdify(x, sol_y, 'numpy')
y_arr = y_func(x_arr)

plt.plot(x_arr, y_arr)
plt.show()

二、复杂隐函数的处理（无解析解） #

你的实际表达式b_sim = -1无法通过符号方法得到解析解，因此需要用数值求解的方式：

1. 固定一个变量的取值范围（比如nu从0.225到1.5）
2. 对每个nu值，用数值方法求解对应的h值
3. 收集所有(nu, h)对后绘图

示例代码：

import sympy as sym
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import root_scalar

# 定义符号变量和表达式
h, nu = sym.symbols('h nu') 
b_sim = sym.cos(sym.pi*sym.sqrt(1 - h)/(2*nu))*sym.cos(sym.pi*sym.sqrt(h + 1)/(2*nu)) - sym.sin(sym.pi*sym.sqrt(1 - h)/(2*nu))*sym.sin(sym.pi*sym.sqrt(h + 1)/(2*nu))/sym.sqrt(1 - h**2)
# 定义方程：b_sim + 1 = 0
eq = b_sim + 1

# 转换为数值函数：输入h和nu，返回方程值
eq_func = sym.lambdify([h, nu], eq, 'numpy')

# 定义针对每个nu值求解h的函数
def solve_h_for_nu(nu_val):
    # 固定nu，生成单变量求解函数
    def func(h_val):
        return eq_func(h_val, nu_val)
    # 尝试在指定区间内找根（brentq方法要求区间端点函数值符号相反）
    try:
        res = root_scalar(func, bracket=[-1.0, 1.0], method='brentq')
        return res.root if res.converged else np.nan
    except:
        return np.nan

# 生成nu的数值数组
nu_arr = np.linspace(0.225, 1.5, 200)
# 求解对应的h值
h_arr = [solve_h_for_nu(nu_val) for nu_val in nu_arr]

# 绘图
plt.plot(nu_arr, h_arr)
plt.xlabel('nu')
plt.ylabel('h')
plt.ylim(-1.1, 1.1)
plt.show()

注意事项 #

• 数值求解时需要根据方程特性调整求解区间和方法，确保函数在区间端点符号相反，才能让brentq这类方法生效
• 若某些nu值无有效根，返回NaN可让Matplotlib自动跳过这些点，避免绘图错误
• 增加nu数组的采样点数量，可让绘图线条更平滑

内容的提问来源于stack exchange，提问作者J. Serra