如何用SymPy求解与theta无关的c1、c2系数？

我有方程 0 = c1*f(x,y,z) + c2*g(x,y,z)，希望用SymPy求解c1、c2作为y,z的函数，已知存在不含theta的解，但现有代码得到的解包含theta。

原尝试代码 #

from sympy import *
# definitions
theta, phi, psi = symbols('theta phi psi')
aa, ab, ac, bb, bc, cc, c1, c2 = symbols('aa ab ac bb bc cc c1 c2')
a, b, c, R = symbols('a b c R')
lab = ['x', 'y', 'z']

# direction cosine matrix
R1 = Matrix([[cos(psi),-sin(psi),0],[sin(psi),cos(psi),0],[0,0,1]])
R2 = Matrix([[cos(theta),0,sin(theta)],[0,1,0],[-sin(theta),0,cos(theta)]])
R3 = Matrix([[cos(phi),-sin(phi),0],[sin(phi),cos(phi),0],[0,0,1]])

DCM = R1*R2*R3
alpha = Matrix([[aa,ab ,ac ], [ab, bb, bc], [ac, bc, cc]])
mu = Matrix([a, b, c])

muLab = DCM*mu
alphaLab = DCM*alpha*DCM.T
#R = aa/cc
# lab frame
chi2 = MutableDenseNDimArray(zeros(27), (3,3,3))
for i in [0, 1, 2]:
    for j in [0, 1, 2]:
        for k in [0, 1, 2]:
            chi2[i,j,k] = alphaLab[i,j]*muLab[k]

norm = 4*pi**2
xxz = simplify(integrate(chi2[0,0,2], (phi, 0, 2*pi), (psi, 0, 2*pi))/norm)
xzx = simplify(integrate(chi2[0,2,0], (phi, 0, 2*pi), (psi, 0, 2*pi))/norm)
zzz = simplify(integrate(chi2[2,2,2], (phi, 0, 2*pi), (psi, 0, 2*pi))/norm)

print(solve((c1*xxz-c2*xzx-zzz),c1,c2))

这段代码无法得到不含theta的解。

简化版问题 #

xxz2 = cc*c*(1/2)*((1+R)*cos(theta)-(1-R)*(cos(theta)**3))
xzx2 = cc*c*(1/2)*(1-R)*(cos(theta)-cos(theta)**3)
zzz2 = cc*c*(R*cos(theta)+(1-R)*(cos(theta)**3))

print(solve((c1*xxz2-c2*xzx2-zzz2),c1,c2,check=False))

同样无法得到不含theta的c1、c2解，但手动指定c1=1/R后：

print(solve(((1/R)*xxz2-c2*xzx2-zzz2), c2))

能正确得到解：

{c2: 1/R + 2}

原因分析 #

• 你只有一个方程，但要解两个变量c1、c2，SymPy默认会将其中一个变量视为自由参数，用另一个变量和theta表示通解，这是线性方程组的常规解形式。
• 你需要的是不含theta的特解，但SymPy的solve函数默认优先输出通解，不会主动筛选满足“解与theta无关”的额外约束的特解。

解决方法 #

方法1：添加“解与theta无关”的约束 #

通过要求c1、c2对theta的导数为0，明确约束解不含theta：

from sympy import *
theta, R, cc, c, c1, c2 = symbols('theta R cc c c1 c2')

xxz2 = cc*c*(1/2)*((1+R)*cos(theta)-(1-R)*(cos(theta)**3))
xzx2 = cc*c*(1/2)*(1-R)*(cos(theta)-cos(theta)**3)
zzz2 = cc*c*(R*cos(theta)+(1-R)*(cos(theta)**3))

# 原方程 + c1、c2与theta无关的约束（导数为0）
eq1 = c1*xxz2 - c2*xzx2 - zzz2
eq2 = diff(c1, theta)
eq3 = diff(c2, theta)

sol = solve((eq1, eq2, eq3), c1, c2, dict=True)
print(sol)

运行结果：[{c1: 1/R, c2: 1/R + 2}]

方法2：提取多项式系数令其为0 #

将方程整理为cos(theta)的多项式，令cos^3(theta)和cos(theta)的系数分别为0，构建方程组求解：

from sympy import *
theta, R, cc, c, c1, c2 = symbols('theta R cc c c1 c2')

xxz2 = cc*c*(1/2)*((1+R)*cos(theta)-(1-R)*(cos(theta)**3))
xzx2 = cc*c*(1/2)*(1-R)*(cos(theta)-cos(theta)**3)
zzz2 = cc*c*(R*cos(theta)+(1-R)*(cos(theta)**3))

eq = c1*xxz2 - c2*xzx2 - zzz2
# 提取cos^3(theta)和cos(theta)的系数
coeff_cos3 = eq.coeff(cos(theta)**3)
coeff_cos = eq.coeff(cos(theta))

sol = solve((coeff_cos3, coeff_cos), c1, c2, dict=True)
print(sol)

运行结果：[{c1: 1/R, c2: 1/R + 2}]

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