sympy 没有给我该符号的确切答案

Question

我有多个方程，我想用

sympy

得到 NT 作为数字，但它没有按我的预期工作。

以下是方程式：

这是我的代码

with open("t.txt","r") as file:
    lines = file.readlines()

t10 = t32 = 0.1
t21 = 3.2
g0=3
g1 = 1
g2 = 1
g3 = 2
R1 = R2 = 0.995
L = 120
l = 5
a_cm = 4059
Eresult32_eV = 4.71 # E3 - E2 eV
Eresult32_nm = 270 # E3 - E2 nm
Eresult10_eV = (4.22-4.78) # E3 - E2 eV
Eresult10_nm = (298-259) # E3 - E2 
K = (1.38*(10**-23))
sigma_mo = (0.92*(10**-19))
p_in = 3.73
a_persent = 1.7
hv_L = ((10**-17)*0.026)
T_R = (10**-6)*0.0544
P_0 = 5
pi = 3.14159
c= 3*(10**8)

def Wpx(z):
    return (75.64*(10**-6)*((1+311.210665*(z-0.0025)**2)**sp.Rational(1,2)))
def Wpy(z):
    return (77.41*(10**-6)*((1+304.239992*(z-0.0025)**2)**sp.Rational(1,2)))
def l_pump(x,y,z):
    # return (2*P_0)/(pi*Wpx(z)*Wpy(z))*sp.exp((((-2*x**2)/(Wpx(z)**2))-((-2*y**2)/(Wpy(z)**2))-a_p*z))
    return ((2*P_0)/(pi*Wpx(z)*Wpy(z)))*sp.exp(((-2*x**2)/(Wpx(z)**2))-((-2*y**2)/(Wpy(z)**2))-(a_p*z))
def Wp (x,y,z):
    return (sigma_mo*l_pump(x,y,z))/hw
    
a_p = 4.043
landa = 455
hw = ((10**-17)*0.0259)

N0,N1,N2,N3,phi,NT = sp.symbols("N0,N1,N2,N3,phi,NT",real=True, positive=True)
NT_array = np.zeros(len(lines),dtype=float)
z_array = np.zeros(len(lines),dtype=float)
x_array = np.zeros(len(lines),dtype=float)
y_array = np.zeros(len(lines),dtype=float)
temp_array=np.zeros(len(lines),dtype=float)

for index,line in enumerate(lines) :
    x,y,z,temp = line.split()
    x,y,z,temp = float(x),float(y),float(z),float(temp)
    z_array[index] = z
    x_array[index] = x
    y_array[index] = y
    temp_array[index] = temp
    
    dtN0 = sp.Eq((-Wp(x,y,z)*N0)+(N1/t10),0)
    dtN1 = sp.Eq((N2-(g2/g1)*N1)*sigma_mo*phi-(N1/t10)+(N2/t21),0)
    dtN2 = sp.Eq(-(N2-(g2/g1)*N1)*sigma_mo*phi+(N3/t32)-(N2/t21),0)
    dtN3 = sp.Eq((Wp(x,y,z)*N0)-(N3/t32),0)
    eq5 = sp.Eq( N0 + N1 + N2 + N3,NT)  # Total population equation
    resulte = sp.solve((dtN0,dtN1,dtN2,dtN3,eq5),(NT))
    print(resulte)

这是我得到的输出

{NT：N0 + N1 + N2 + N3}

{NT：N0 + N1 + N2 + N3}

{NT：N0 + N1 + N2 + N3}

。 。 .

所以我期望

NT

是一个数字，而不是一个方程。我还尝试更改代码以获取

N0,N1,N2,N3

并使用

sum()

来获取

NT

作为公式

N0+N1+N2+N3=NT

，但随后我得到了此输出

{N0: 1.02832264875651e-220*NT, N1: 0.0294117647058824*NT, N2: 0.941176470588235*NT, N3: 0.0294117647058824*NT}
{N0: 6.58641650153513e-226*NT, N1: 0.0294117647058824*NT, N2: 0.941176470588235*NT, N3: 0.0294117647058824*NT}
{N0: 8.37786843284486e-219*NT, N1: 0.0294117647058824*NT, N2: 0.941176470588235*NT, N3: 0.0294117647058824*NT}

那有什么问题呢。我怎样才能得到NT作为号码？

Answer 1

让我们包括您的

dphi

方程，因为给定

NT

，这将确定

Ni

其他值之间的关系。让我们象征性地这样做：

from sympy import *
N1, N3, phi, N0, t32, g21, c, l, alpha, NT, t10, sigma, L, N2, sigma_mo, wp, R12, t21 = symbols('N1, N3, phi, N0, t32, g21, c, l, alpha, NT, t10, sigma, L, N2, sigma_mo, w_p, R12, t21')
dtN0 = Eq((-wp*N0)+(N1/t10),0)
dtN1 = Eq((N2-g21*N1)*sigma_mo*phi-(N1/t10)+(N2/t21),0)
dtN2 = Eq(-(N2-g21*N1)*sigma_mo*phi+(N3/t32)-(N2/t21),0)
dtN3 = Eq((wp*N0)-(N3/t32),0)
eq5 = Eq( N0 + N1 + N2 + N3,NT)  # Total population equation
dphi = Eq((N2-g21*N1)*sigma*l*c*phi/L-phi*c*(1-R12)/2/L-2*alpha*l*phi,0)
seqs = [(i.lhs - i.rhs).as_numer_denom()[0] for i in [dtN0,dtN1,dtN2,dtN3,dphi,eq5]]
sol = solve(seqs, (N0,N1,N2,N3,phi), dict=True)

有两种解决方案：一种是

phi=0

，另一种是非零值。它们都满足方程组。

>>> [[i.subs(s).cancel() for i in seqs] for s in sol]
[[0,0,0,0,0,0], [0,0,0,0,0,0]]

以下是每种情况下必须指定的符号：

>>> [list(ordered(Dict(s).free_symbols)) for s in sol]
[[N0, N1, N2, N3, NT, phi, t10, t21, t32, w_p], [L, N0, N1, N2, N3, NT, R12, alpha, c, g21, l, phi, sigma, sigma_mo, t10, t21, t32, w_p]]

在

phi=0

的情况下，所有

Ni

都与

NT

成比例，因此您可以让

NT

为 1，然后

Ni

值代表实际

NT

值的分数。在

Ni

取决于

NT

但不成比例的非零情况下，情况并非如此。如果您打印解决方案，您将看到这个，例如

>>> sol[1][N0]
(-4*L*alpha*l + 2*NT*c*l*sigma + R12*c - c)/(2*c*l*sigma*(g21*t10*w_p + t10*w_p + t32*w_p + 1))

我们可以通过寻找解中不能为 0 的东西（即分母）来检查常数方面的约束：

from sympy.solvers.solvers import denoms
>>> [denoms(Dict(s)) for s in sol]
[{t10*w_p + t21*w_p + t32*w_p + 1}, {-4*L*alpha*l + R12*c - c, 2, c, l, sigma, g21*t10*w_p + t10*w_p + t32*w_p + 1, sigma_mo, t21}]

对于每个解决方案，这些值不能为 0。我鼓励您使用分母来测试您的常量：如果对于给定值，它们都不为 0，那么您将有一个解决方案（看起来）。 （您还可以包含更多定义

w_p

的符号；我也不清楚

sigma

和

sigma_mo

之间的关系。）