如何用SymPy构建x与y的无穷级数函数并修复代码报错？

问题根源 #

原代码报错function object cannot be interpreted as integer的核心原因：

• phi_exact中调用total(x,y,N)时N未定义，且试图将Python函数传入SymPy极限计算接口，但SymPy仅支持符号表达式的极限运算
• 混用SymPy符号计算与NumPy数值函数（如np.sin搭配sym.pi），导致类型冲突
• 嵌套函数不符合导师要求，且逻辑混乱

修复方案与代码 #

方案1：符号计算求无穷级数精确解 #

利用SymPy的sym.summation直接计算无穷级数，避免手动循环和错误的极限调用：

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

# 物理常数
e_0 = 8.85e-12
k_e = 1/(4*sym.pi*e_0)
pi = sym.pi

# 变分法近似解（数值函数）
A = 5/(4*e_0)
def phi_var(x, y):
    return A * x * (1 - x) * y * (1 - y)

# 精确解：用SymPy定义无穷级数符号表达式
x_sym, y_sym, m_sym = sym.symbols('x y m', integer=True, positive=True)
n_sym = 2*m_sym + 1  # 仅取奇数项
term = (sym.sin(n_sym * pi * x_sym) / n_sym**3) * (1 - sym.cosh(n_sym * pi * (y_sym - 1/2)) / sym.cosh(n_sym * pi / 2))
exact_series = sym.summation(term, (m_sym, 1, sym.oo))
phi_exact_sym = (16 * k_e / pi**2) * exact_series

# 转成数值函数，方便绘图计算
phi_exact_num = sym.lambdify((x_sym, y_sym), phi_exact_sym, 'numpy')

# 测试计算
print("变分解(1,0):", phi_var(1, 0))
print("精确解(1,0):", phi_exact_num(1, 0))

方案2：大N数值近似无穷级数 #

如果不需要严格符号解，利用级数收敛快（分母为n³）的特点，取足够大的N（如100）近似：

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

e_0 = 8.85e-12
k_e = 1/(4*np.pi*e_0)  # 统一用NumPi做数值计算
pi = np.pi

# 变分法近似解
A = 5/(4*e_0)
def phi_var(x, y):
    return A * x * (1 - x) * y * (1 - y)

# 精确解：大N近似无穷级数
def phi_exact_num(x, y, N=100):
    coeff = 16 * k_e / pi**2
    total = 0.0
    for m in range(1, N+1):
        n = 2*m + 1
        trig = np.sin(n * pi * x)
        denom = n**3
        cosh_numer = np.cosh(n * pi * (y - 0.5))
        cosh_denom = np.cosh(n * pi * 0.5)
        ratio = cosh_numer / cosh_denom
        total += (trig / denom) * (1 - ratio)
    return coeff * total

# 测试计算
print("变分解(1,0):", phi_var(1, 0))
print("精确解(1,0):", phi_exact_num(1, 0))

绘图对比示例 #

生成网格点计算两个解的值，绘制热力图对比：

# 生成网格
x = np.linspace(0, 1, 50)
y = np.linspace(0, 1, 50)
X, Y = np.meshgrid(x, y)

# 计算解的值
var_sol = phi_var(X, Y)
exact_sol = phi_exact_num(X, Y)  # 对应方案2；若用方案1则直接调用phi_exact_num(X,Y)

# 绘图
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 6))

# 变分解热力图
im1 = ax1.imshow(var_sol, extent=[0,1,0,1], origin='lower', cmap='viridis')
ax1.set_title('变分法近似解')
plt.colorbar(im1, ax=ax1)

# 精确解热力图
im2 = ax2.imshow(exact_sol, extent=[0,1,0,1], origin='lower', cmap='viridis')
ax2.set_title('精确解（无穷级数近似）')
plt.colorbar(im2, ax=ax2)

# 差值热力图
diff = np.abs(var_sol - exact_sol)
im3 = ax3.imshow(diff, extent=[0,1,0,1], origin='lower', cmap='Reds')
ax3.set_title('近似解与精确解差值')
plt.colorbar(im3, ax=ax3)

plt.tight_layout()
plt.show()

关键修改说明 #

1. 符号计算版本：用sym.summation直接定义无穷级数，避免手动循环和错误的极限调用，同时将符号表达式转为数值函数用于绘图
2. 数值近似版本：统一用NumPy做数值计算，去掉嵌套函数，取N=100足够近似无穷级数（因n³衰减，前100项已足够精确）
3. 移除原代码中未定义的变量引用和函数嵌套，符合导师要求

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