如何在函数中忽略符号变量a，对其系数执行数值运算？

问题描述 #

编写了如下Python代码，试图用sympy符号变量a计算晶格势能，但因为符号变量无法执行round、sort等数值操作，报错Cannot round symbolic expression和cannot determine truth value of Relational：

import numpy as np
import scipy as sp
import sympy as smp
import matplotlib.pyplot as plt
from scipy.misc import derivative
from sympy import Eq,solve
from itertools import groupby 

absulom= 2.844 #[mev]
sigma= 0.34 #[nm]

def potential(a):
    a1= np.array([3/2*a,3.**0.5/2.*a])
    a2= np.array([3/2*a,-3.**0.5/2.*a])
    nmax = 10
    coordsA =[ i * a1 + j* a2 for i in range(-nmax,nmax) for j in range(-nmax,nmax)]
    coordsB= [i * a1 + j* a2 + [a,0.] for i in range(-nmax,nmax) for j in range(-nmax,nmax)]
    distance1=[((- coordsA[k][0])**2+(-coordsA[k][1])**2)**0.5 for k in range(len(coordsA)-1)]
    distance1= [round(num, 5.) for num in distance1]
    distance2=[((- coordsB[s][0])**2+(-coordsB[s][1])**2)**0.5 for s in range(len(coordsB)-1)]
    distance2= [round(num, 5.) for num in distance2]
    distance=distance1+distance2
    distance.sort()
    distance.pop(0)
    grouped_neighbors=[list(j) for i, j in groupby(distance)]
    woutrep=[grouped_neighbors[v][0]for v in range(31)]
    U=0
    for t in range(30):
        V=2*absulom*((sigma/grouped_neighbors[t][0])**12-(sigma/grouped_neighbors[t][0])**6)*np.size(grouped_neighbors[t])
        U=U+V

a = smp.symbols('a', real=True)
f = potential(a)

解决方案 #

核心思路是：所有晶格点到原点的距离都是a乘以一个无量纲系数，因此可以预先用数值计算（把a=1）得到这些系数，完成排序、去重、分组操作，再把结果和符号a结合进行势能计算，避开对符号表达式做数值操作。

修改后的代码 #

import numpy as np
import sympy as smp
from itertools import groupby 

absulom= 2.844  #[mev]
sigma= 0.34     #[nm]

# 预先计算距离系数和对应邻居数量（把a=1处理）
def precompute_coeffs():
    a = 1.0  # 用数值1代替符号a，计算系数
    a1 = np.array([3/2*a, 3.**0.5/2.*a])
    a2 = np.array([3/2*a, -3.**0.5/2.*a])
    nmax = 10
    
    coordsA = [i * a1 + j * a2 for i in range(-nmax, nmax) for j in range(-nmax, nmax)]
    coordsB = [i * a1 + j * a2 + [a, 0.] for i in range(-nmax, nmax) for j in range(-nmax, nmax)]
    
    # 计算距离系数（此时a=1，结果就是系数）
    distance1 = [((-coord[0])**2 + (-coord[1])**2)**0.5 for coord in coordsA[1:]]  # 跳过原点
    distance2 = [((-coord[0])**2 + (-coord[1])**2)**0.5 for coord in coordsB]
    
    distance = distance1 + distance2
    distance = [round(num, 5) for num in distance]
    distance.sort()
    
    # 分组并提取每组的系数和数量
    grouped = [list(g) for _, g in groupby(distance)]
    # 取前30组（对应原代码的range(30)）
    coeff_counts = [(g[0], len(g)) for g in grouped[:30]]
    return coeff_counts

# 预先计算好系数和数量
coeff_counts = precompute_coeffs()

def potential(a):
    U = 0
    for coeff, count in coeff_counts:
        # 实际距离是 coeff * a，代入势能公式
        r = coeff * a
        V = 2 * absulom * ((sigma / r)**12 - (sigma / r)**6) * count
        U += V
    return U

# 符号变量计算
a = smp.symbols('a', real=True, positive=True)  # 加positive确保r不为负
f = potential(a)
# 可以化简结果
print(smp.simplify(f))

关键修改说明 #

• 预计算系数：用a=1先算出所有距离的无量纲系数，完成round、sort、groupby这些数值操作，得到每个系数对应的邻居数量。
• 符号计算分离：在势能函数中，直接用预计算的系数乘以符号a得到实际距离，再代入公式计算，全程只对符号a进行代数运算，避免了数值操作和符号操作的冲突。
• 优化细节：给符号a加上positive=True，避免后续运算中出现根号下负数的问题；跳过原点（coordsA[1:]），和原代码的pop(0)逻辑一致。

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