I am new to both python and the sympy module so some of my code may be a bit 'primitive'. This is also my first use of Stackoverflow so forgive any use errors. I have included an entire test function which can be called using TEST() [this is a simplified version of my actual work, but displays the problem I am getting]. Essentially I have three equations:
(1) lamb=h/sqrt(k*T)
(2) lamb=d
(3)d**3=k*T/P
and I wish to amalgamate them and solve for T. This is quite easy to do by hand without using the sympy package, but, as I say, I am trying to learn it. When I run the code, instead of getting the answer T=h**(6/5)P**(2/5)/k, I get a series of complex numbers. I have no problem if, instead of solving for T, I solve for h. Moreover, if I change equation (3) to, say d**4=k*T/P, the result I get is correct. The problem seems to arise out of the fact that in trying to solve the equation, the powers of k and T (in the equations) are the same. In fact, if I replace k*T by the single variable kT throughout, I get a sensible solution.
Can anyone shed some light on what I am doing wrong and what I should be doing instead?
MORE INFO: In fact the correct solution is present! It's the [4]th answer. But why am I getting the other (complex number) answers (and how do I select the real valued one)?
MORE THOUGHTS (16/-6/20 18:15): Presumably it is because k and T independently can have 5 roots (because of the 6/5 power). So is there a way of saying only real roots should be considered?
PaulJ
def TEST():
import sympy
from sympy import symbols,Symbol,root,integrate,diff,Derivative,N,I,var,latex
from sympy import sympify,Function,Eq,solve
from sympy import pi,exp
print("")
print("TEST")
#Declare other variables (avoid C,O,S,I,N,E, and Q)
var('d h k T P lamb kT',positive=True)
Eq1=Eq(lamb,h/(root(T,2,0)*root(k,2,0)))
print("Eq1=",Eq1)
Eq2=Eq(lamb,d)
print("Eq2=",Eq2)
Eq3=Eq(d**3,k*T/P)
print("Eq3=",Eq3)
res_d=solve(Eq3,d)[0]
print("res_d=",res_d)
Eq4=Eq2.subs(d,res_d)
print("Eq4=",Eq4)
res_lamb=solve(Eq4,lamb)[0]
print("res_lamb=",res_lamb)
Eq5=Eq1.subs(lamb,res_lamb)
print("Eq5=",Eq5)
res_P=solve(Eq5,P)[0]
print("res_P=",res_P)
Eq6=Eq(res_P,P)
print("Eq6=",Eq6)
print("Try solving for (this works) h=",solve(Eq6,h)[0])
print("Try solving for (this fails miserably) T=",solve(Eq6,T)[0])
print("and in fact, the full miserable solution is...")
res_T=solve(Eq6,T)[0]
print("res_T=",solve(Eq6,T))
print("")
print("The problem seems to be that the powers of T and k in Eq6 are the same.")
print("If I change the Eq3 so that it's d**4 = k*T/P I get a sensible solution")
print("")
print("Now try with k*T replaced by kT throughout")
Eq1a=Eq(lamb,h/(root(kT,2,0)))
print("Eq1a=",Eq1a)
Eq2a=Eq(lamb,d)
print("Eq2a=",Eq2a)
Eq3a=Eq(d**3,kT/P)
print("Eq3a=",Eq3a)
resa_d=solve(Eq3a,d)[0]
print("resa_d=",resa_d)
Eq4a=Eq2.subs(d,resa_d)
print("Eq4a=",Eq4a)
resa_lamb=solve(Eq4a,lamb)[0]
print("res_lamb=",resa_lamb)
Eq5a=Eq1a.subs(lamb,resa_lamb)
print("Eq5a=",Eq5a)
resa_P=solve(Eq5a,P)[0]
print("res_P=",resa_P)
Eq6a=Eq(resa_P,P)
print("Eq6a=",Eq6a)
print("Try solving for (this works) h=",solve(Eq6a,h)[0])
print("Try solving for (this work!) kT=",solve(Eq6a,kT)[0])
return
