如何绘制双变量SymPy函数？C_EF函数实现后绘图不符问题求助

核心问题分析 #

从代码逻辑和预期/实际绘图的差异来看，以下几个关键点可能导致结果不符：

1. 符号与数值混合运算的冲突 #

用SymPy定义了符号变量wd，但dlt直接通过.evalf()转为数值常量，后续C_EF表达式同时混合符号变量和数值，可能导致SymPy解析绘图逻辑时出现异常，无法正确反映w/δ与C_EF的函数关系。

2. 求和收敛判断的错误 #

求和终止条件中使用math.fabs(term_1)，但term_1是SymPy表达式对象，math.fabs无法正确获取其数值绝对值，可能导致求和提前终止或无法收敛到足够精度，最终影响C_EF的计算结果。

3. 复数运算的精度问题 #

计算beta_n时的复数平方根处理、实部虚部的提取逻辑，以及abs(beta_n)**6的计算方式，可能与原公式的数学定义存在偏差，SymPy对复数的符号运算精度有限，容易引入误差。

4. 变量作用域与逻辑一致性 #

l是全局定义的3*w_magnet，但函数P_calc接收参数w，计算lmd_n时使用的是参数w，这种全局+局部变量的混合方式容易引发逻辑混乱。

修正后的代码示例 #

### IMPORTS AND ABBREVATIONS ###
import numpy as np
import scipy
from scipy.special import sinh, cosh, sin, cos
import matplotlib.pyplot as plt

mu_r = 1.05  # 相对磁导率
mu_0 = 4 * scipy.pi * 10 ** -7  # 真空磁导率
mu = mu_0 * mu_r  # 绝对磁导率 (H/m)
sigma = 7.1 * 10 ** 5  # 电导率 (S/m，NdFeB磁铁)
f = 2500  # 载波频率 (Hz)
B_h_s = 1.5  # 峰值磁通密度 (T)

dAG = 2.5 * 10 ** -3  
w_magnet = 60 * 10 ** -3  
h = 5 * 10 ** -3  
l = 3 * w_magnet  
pi = scipy.pi

def P_calc(w):
    # 计算dlt数值
    dlt = np.sqrt(1 / (pi * mu * sigma * f)) * np.sqrt((h + dAG) / h)

    def compute_converged_sum():
        tolerance = 1e-6
        max_iterations = 1000
        sum_total = 0.0
        n = 0
        while True:
            lmd_n = (2 * n + 1) * pi / w
            # 复数beta_n计算
            beta_sq = lmd_n ** 2 + 2 * 1j / (dlt ** 2)
            beta_n = np.sqrt(beta_sq)
            re_beta = beta_n.real
            im_beta = beta_n.imag
            
            denominator = (2 * n + 1) ** 5 * (abs(beta_n) ** 6) * (cosh(re_beta * l) + cos(im_beta * l))
            # 计算两项并累加
            term1 = ((lmd_n ** 2 - 2 * im_beta ** 2) * re_beta * lmd_n ** 3 * sinh(re_beta * l)) / denominator
            term2 = ((lmd_n ** 2 + 2 * re_beta ** 2) * im_beta * lmd_n ** 3 * sin(im_beta * l)) / denominator
            sum_total += term1 + term2
            
            # 双重收敛判断，确保两项都足够小
            if abs(term1) < tolerance and abs(term2) < tolerance or n >= max_iterations:
                break
            n += 1
        return sum_total.real  # 理论上求和结果应为实数

    sum_C_EF = compute_converged_sum()
    
    # 生成w/δ的数值范围，避开0点避免除零错误
    wd_vals = np.linspace(0, 6, 1000)
    wd_safe = np.where(wd_vals == 0, 1e-12, wd_vals)
    
    # 计算C_EF各部分
    C_EF_sinh_cosh = (cosh(wd_safe) + cos(wd_safe)) / (sinh(wd_safe) - sin(wd_safe))
    C_EF_term = ((32 * w) / ((pi ** 5) * l)) * (wd_safe ** 3)
    C_EF = 1 - C_EF_term * C_EF_sinh_cosh * sum_C_EF

    # 绘图
    plt.plot(wd_vals, C_EF)
    plt.ylabel('C_EF')
    plt.xlabel('w/δ')
    plt.grid(True)
    plt.show()

    # 补充P_e的计算逻辑（根据原公式完善）
    # P_e = ...
    return None

P_calc(w_magnet)

修正要点说明 #

• 改用numpy/scipy进行全数值计算，避免SymPy符号运算的解析误差，绘图逻辑更稳定。
• 合并求和循环，统一收敛判断，使用复数绝对值确保求和精度。
• 处理wd=0时的除零问题，避免程序报错。
• 提取求和结果的实部，符合C_EF为实数的物理意义。

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