As Sage fails to solve the non-linear (2nd order polynomial) systems am I
interested in, I tried to force the substitution of variables in the system
to obtain an answer.
In the example given below, the system corresponds to :
g1 = 0
g2 = 0
where the unknowns are P1m and P2p. The others variables, P1p, M1, M2 and
gm (s=0) are considered as parameters.
g1 and g2 being polynomials of P1m and P2p, here is what I do manually :
1° I compute the P2p (named P2p_root) root of g1. P2p_root thus depends on
P1m.
2° I evaluate g2 with P2p = P2P_root. Now, g2(P2p=P2p_root) only depends on
P1m (and parameters P1p, M1, M2 and gm).
3° I compute the P1m root of g2(P2p=P2P_root).
4° I finally evaluate P2p_root with the value obtained at step 3° for P1m.
Il works well for the first order system. However, for the second order
system the resolution fails during step 3. Indeed, the P1m root of
g2(P2p=P2P_root) contains the variable P1m, which sounds absurd to me.
I don't understand why the command 'solve(g2(P2p=P2p_root)==0,P1m)' returns
a result that depends on P1m !
Can please anybody help me? I really need to get the solutions of this
problem and it is driving me crazy !
Best regards.
#################
## My code
#################
Pp = var('Pp')
Pm = var('Pm')
s = var('s')
M = var('M')
gm = var('gm')
#P1m and P2p are the unknowns
P1p = var('P1p')
P1m = var('P1m')
P2p = var('P2p')
P2m = var('P2m')
M1 = var('M1')
M2 = var('M2')
assume(P1p > 0)
assume(M1 > 0)
assume(M2 > 0)
# First order
# -----------
print "First order resolution"
print "----------------------"
f1 = (1+1/M)*Pp + (1-1/M)*Pm - s
f2 = ( (gm-1)*(1+M)*Pp + (gm-1)*(1-M)*Pm + s ) / (1+(gm-1)/2*M^2)
g1 = f1(Pp=P1p,Pm=P1m,s=0,M=M1) - f1(Pp=P2p,Pm=0,s=0,M=M2)
g2 = f2(Pp=P1p,Pm=P1m,s=0,M=M1) - f2(Pp=P2p,Pm=0,s=0,M=M2)
P2p_root = solve(g1==0,P2p)[0].rhs()
P1m_root = solve(g2(P2p=P2p_root)==0,P1m)[0].rhs()
P2p_root = P2p_root(P1m=P1m_root)
print "P1m = ", P1m_root.factor()
print "P2p = ", P2p_root.factor()
# Second order
# ------------
print
print "Second order resolution"
print "-----------------------"
f1 = (1+1/M)*Pp + (1-1/M)*Pm - s + s^2/2 -
((gm-1)*M-2)/(2*M)*Pp^2 \
- ((gm-1)*M+2)/(2*M)*Pm^2 - (gm-1)*Pm*Pp - (1+1/M)*Pp*s - (1-1/M)*Pm*s
f2 = ( (gm-1)*(1+M)*Pp + (gm-1)*(1-M)*Pm + s + s^2/2 -
2*(gm-1)*Pp*Pm \
+ (gm-1)*Pp*s + (gm-1)*Pm*s ) / (1+(gm-1)/2*M^2)
g1 = f1(Pp=P1p,Pm=P1m,s=0,M=M1) - f1(Pp=P2p,Pm=0,s=0,M=M2)
g2 = f2(Pp=P1p,Pm=P1m,s=0,M=M1) - f2(Pp=P2p,Pm=0,s=0,M=M2)
P2p_root = solve(g1==0,P2p)[0].rhs()
P1m_root = solve(g2(P2p=P2p_root)==0,P1m)[0].rhs()
P2p_root = P2p_root(P1m=P1m_root)
print "P1m = ", P1m_root.factor()
print "P2p = ", P2p_root.factor()
--
To post to this group, send email to sage-support@googlegroups.com
To unsubscribe from this group, send email to
sage-support+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/sage-support
URL: http://www.sagemath.org
