Solve for variable but variable is still in answer

I try to solve for variable P3. Sage gives me an answer for P3 but it still contains P3 on the right hand side (see last equation after the solve in the beginning of sqrt?

sage: A1, A2, A3, P1, P2, P3, u1, u2, u3, r1, r2, r3, g, C, M1, M2, M3 = var('A1 A2 A3 P1 P2 P3 u1 u2 u3 r1 r2 r3 g C M1 M2 M3')
sage: eq1=A3*(P1-P3)==A3*r3*u3^2-(A2*r2*u2^2+A1*r1*u1^2)
sage: eq2=A3*u3*r3==A2*u2*r2+A1*r1*u1
sage: eq5a=solve([eq2],r3)
sage: eq5b=u3^2+2*g/(g-1)*(P3/(eq5a[0].rhs()))-2*C==0
sage: eq6=solve([eq5b],u3)
sage: eq7a=eq1/A3
sage: eq7b=eq7a.subs(-A3*r3*u3^2 == -(A2*u2*r2+A1*r1*u1)*u3)
sage: eq7b
P1 - P3 == -(A1*r1*u1^2 + A2*r2*u2^2 + (A3*P3*g + sqrt(A3^2*P3^2*g^2 + 2*((g^2 - 2*g + 1)*A1^2*r1^2*u1^2 + 2*(g^2 - 2*g + 1)*A1*A2*r1*r2*u1*u2 + (g^2 - 2*g + 1)*A2^2*r2^2*u2^2)*C))*(A1*r1*u1 + A2*r2*u2)/(A1*(g - 1)*r1*u1 + A2*(g - 1)*r2*u2))/A3

sage: eq7c=eq7b.subs(u3==eq6[0].rhs())
sage: eq7=solve([eq7c],P3)
sage: eq7
[P3 == (A1*r1*u1^2 + A2*r2*u2^2 + A3*P1 - (A1*r1*u1^2 + A2*r2*u2^2 + A3*P1)*g - sqrt(A3^2*P3^2*g^2 + 2*((g^2 - 2*g + 1)*A1^2*r1^2*u1^2 + 2*(g^2 - 2*g + 1)*A1*A2*r1*r2*u1*u2 + (g^2 - 2*g + 1)*A2^2*r2^2*u2^2)*C))/A3]

edit retag close merge delete

It looks like it sees P3 and P3 as different variables?

( 2014-08-28 22:39:48 +0200 )edit

Trying to recreate this gave me an error because eq2 was not defined.

( 2014-08-29 07:47:53 +0200 )edit

I forgot to copy that line or something, sorry about that

( 2014-08-29 15:48:42 +0200 )edit

Sort by » oldest newest most voted

First, you should notice that eq6 is a list of two solutions for u3:

sage: len(eq6)
2
sage: eq6
[u3 == -(A3*P3*g + sqrt(A3^2*P3^2*g^2 + 2*((g^2 - 2*g + 1)*A1^2*r1^2*u1^2 + 2*(g^2 - 2*g + 1)*A1*A2*r1*r2*u1*u2 + (g^2 - 2*g + 1)*A2^2*r2^2*u2^2)*C))/(A1*(g - 1)*r1*u1 + A2*(g - 1)*r2*u2), u3 == -(A3*P3*g - sqrt(A3^2*P3^2*g^2 + 2*((g^2 - 2*g + 1)*A1^2*r1^2*u1^2 + 2*(g^2 - 2*g + 1)*A1*A2*r1*r2*u1*u2 + (g^2 - 2*g + 1)*A2^2*r2^2*u2^2)*C))/(A1*(g - 1)*r1*u1 + A2*(g - 1)*r2*u2)]


You can choose the first or the second as follows:

sage: eq6[0]
u3 == -(A3*P3*g + sqrt(A3^2*P3^2*g^2 + 2*((g^2 - 2*g + 1)*A1^2*r1^2*u1^2 + 2*(g^2 - 2*g + 1)*A1*A2*r1*r2*u1*u2 + (g^2 - 2*g + 1)*A2^2*r2^2*u2^2)*C))/(A1*(g - 1)*r1*u1 + A2*(g - 1)*r2*u2)
sage: eq6[1]
u3 == -(A3*P3*g - sqrt(A3^2*P3^2*g^2 + 2*((g^2 - 2*g + 1)*A1^2*r1^2*u1^2 + 2*(g^2 - 2*g + 1)*A1*A2*r1*r2*u1*u2 + (g^2 - 2*g + 1)*A2^2*r2^2*u2^2)*C))/(A1*(g - 1)*r1*u1 + A2*(g - 1)*r2*u2)
sage: eq6[0] == eq6[1]
False


Second, you found the limitation of the solve function, which is not very powerful, so it is sometimes not able to find a solution of the form P3 = .... Youe equation involves squares roots and solve does not handle such equations very well. A simpler example is:

sage: eq = x == sqrt(x)
sage: eq.solve(x)
[x == sqrt(x)]


EDIT It seems that sympy is better than Maxima in solving this, so you can do:

sage: import sympy
sage: sympy.solve(eq7c.rhs()-eq7c.lhs(), P3)
[(A3*(A1*r1*u1**2 + A2*r2*u2**2 + A3*P1) - sqrt(A3**2*(-2*A1**2*C*g**2*r1**2*u1**2 + 2*A1**2*C*r1**2*u1**2 + A1**2*g**2*r1**2*u1**4 - 4*A1*A2*C*g**2*r1*r2*u1*u2 + 4*A1*A2*C*r1*r2*u1*u2 + 2*A1*A2*g**2*r1*r2*u1**2*u2**2 + 2*A1*A3*P1*g**2*r1*u1**2 - 2*A2**2*C*g**2*r2**2*u2**2 + 2*A2**2*C*r2**2*u2**2 + A2**2*g**2*r2**2*u2**4 + 2*A2*A3*P1*g**2*r2*u2**2 + A3**2*P1**2*g**2)))/(A3**2*(g + 1)),
(A3*(A1*r1*u1**2 + A2*r2*u2**2 + A3*P1) + sqrt(A3**2*(-2*A1**2*C*g**2*r1**2 ...
more

And there are no other options to solve this? I know there are two solutions (quadratic rule) but I take only one with the eq6[0].rhs() right?

( 2014-08-29 19:21:31 +0200 )edit

yes, eq6 is a Python list that contains all solutions. There are two of them, to get the second one, you just have to look at eq6[1].rhs()

( 2014-08-30 00:23:13 +0200 )edit

Why do you subtract eq7c rhs from lhs **=^

( 2014-09-02 17:18:30 +0200 )edit

@wzawzdb, from the Sympy solvers documentation Use solve() to solve algebraic equations. We suppose all equations are equaled to 0, so solving x**2 == 1 translates into the following code: Run code block in SymPy Live

>>> from sympy.solvers import solve
>>> from sympy import Symbol
>>> x = Symbol('x')
>>> solve(x**2 - 1, x)
[-1, 1]

( 2015-02-25 22:52:22 +0200 )edit