# Solve return an implicit equation

The following solve() resists to any argument (to_poly_solve... calling giac sympy or anything else) in the way that it return an implicit solution in y

var('y')
A=n((1/2)*((1/4)*log(35000*0.25+100)+(3/4)*log(100)))
B=n((1/2)*((1/2)*log(35000*0.25+100)+(1/2)*log(100)))
C=(1/2)*A+ (1/2)*B
D=(1/2)*((1/2)*log(35000*0.25+100-y)+(1/2)*log(100-y))

sol=solve(D-C==0,y)[0]
sol


Of course I can find the solution acting the following way, but I find it after a long search and it seems inelegant.

sol1=e^sol.rhs()-e^sol.lhs()
sol21=n(solve(sol1==0,y)[0].rhs())
sol22=n(solve(sol1==0,y)[1].rhs())
show(sol21)
show(sol22)


Is there an other way ?

edit retag close merge delete

Sage is not smart enough to realize that this equation reduces to a quadratic one. If you don't want to help it, the best shot would be solving the equation numerically via find_root.

( 2021-12-06 18:29:48 +0200 )edit

Sort by » oldest newest most voted

EDIT : The solution below is dubious. See the latter addition...

Not obvious, but there is a way to do this symbolically. Let's get rid of the numerical constants, in order to be able to follow the computations ; we note 35000=K, 0.24=1/4 acd 100=Cent ;

var("K, y")
var("Cent", latex_name="\\mathrm{Cent}")
A=(1/2*(1/4*log(K/4+100)+3/4*log(Cent)))
B=(1/2*(1/2*log(K/4+100)+1/2*log(Cent)))
C=(A+B)/2
D=1/2*(1/2*log(K/4-y)+1/2*log(Cent-y))


Let Eq the equation whose solutions are sought :

sage: Eq=D-C==0 ; Eq
1/4*log(Cent - y) - 5/16*log(Cent) + 1/4*log(1/4*K - y) - 3/16*log(1/4*K + 100) == 0


W get rid of the log fractions (which correspond to fractional powers...) by multiplying my the lcm of their denominators, and push all the "negative" elements of the lhs to the rhs. BEWARE : This multiplication may introduce spurious solutions here !

sage: Eq2=(16*Eq+5*log(Cent)+3*log(K/4+100)).log_simplify() ; Eq2
log(1/256*(Cent - y)^4*(K - 4*y)^4) == log(1/64*Cent^5*(K + 400)^3)


The left and right hands of this equation are equal iif their exponential are equal. Therefore :

sage: Eq3=Eq2.operator()(*map(exp, Eq2.operands())) ; Eq3
1/256*(Cent - y)^4*(K - 4*y)^4 == 1/64*Cent^5*(K + 400)^3


is equivalent to Eq2, and Eq.

Now, we have a vanilla polynomial equation (univariate after reintroduction of the numerical constants), whose roots are trivial to get with Sage :

sage: Sol=(Eq3.lhs()-Eq3.rhs()).subs([Cent==100,K==35000]).roots(ring=QQbar) ; Sol
[(66.7703195987588?, 1),
(8783.22968040125?, 1),
(99.8713723622881? - 33.35634289423430?*I, 1),
(99.8713723622881? + 33.35634289423430?*I, 1),
(133.4869739368500? + 0.?e-137*I, 1),
(8716.51302606315? + 0.?e-68*I, 1),
(8750.12862763772? - 33.35634289423430?*I, 1),
(8750.12862763772? + 33.35634289423430?*I, 1)]


Here, we must filter the spurious solutions we may have introduced before. The best test would be to test the substituted value to zero by is_zero(). This fails for reasons to be explored in a separate ticket. Filtering on a numerical value close to 0 may be "good enough" :

sage: SOL = [u[0] for u in Sol if (D-C).subs({Cent:100, K:35000}).subs(y==u[0]).abs().n()<1e-10] ; SOL
[66.7703195987588?]


BTW :

sage: SOL[0].radical_expression()
-25*sqrt(4*354^(3/4)*sqrt(2) + 29929) + 4425


EDIT 2 : A simpler shortcut leads to doubt the above solution. Working directly on the numerical values :

sage: %cpaste
Pasting code; enter '--' alone on the line to stop or use Ctrl-D.
:var('y')
:A=n((1/2)*((1/4)*log(35000*0.25+100)+(3/4)*log(100)))
:B=n((1/2)*((1/2)*log(35000*0.25+100)+(1/2)*log(100)))
:C=(1/2)*A+ (1/2)*B
:D=(1/2)*((1/2)*log(35000*0.25+100-y)+(1/2)*log(100-y))
:
:--
y
sage: Eq=D==C ; Eq
1/4*log(-y + 8850.00000000000) + 1/4*log(-y + 100) == 3.14314807149665
sage: Eq2 = Eq.operator()(*map(lambda u:u.subs(log(w0)+log(w1)==log(w0*w1)).exp(), (4*Eq).operands())) ; Eq2
(y - 100)*(y - 8850.00000000000) == 288540.947130305


Note that by multiplying the equation, we may have introdiced spurious solutions (see above)...

sage: Sol=Eq2.solve(y) ; Sol
[y == -2/26745*sqrt(3474396435470430) + 4475, y == 2/26745*sqrt(3474396435470430) + 4475]
sage: [u.rhs().n() for u in Sol]
[67.1472407610636, 8882.85275923894]


However, The second solution seems spurious :

sage: [(Eq.lhs()-Eq.rhs()).subs(u).n() for u in Sol]
[4.44089209850063e-16, 5.55111512312578e-16 + 1.57079632679490*I]


I may have erred in my transcription of the original equations. I still think that the principle is sound.

EDIT 3 : A better transcription of the numerical constants gives the expected result. The following code :

# K=35000*0.25
# Data
var('y, K, Cent')
A=(1/2)*((1/4)*log(K+Cent)+(3/4)*log(Cent))
B=(1/2)*((1/2)*log(K+Cent)+(1/2)*log(Cent))
C=(1/2)*A+ (1/2)*B
D=(1/2)*((1/2)*log(K+Cent-y)+(1/2)*log(Cent-y))
# Question
Eq=C==D ; print("Eq = ", Eq)
# Get rid of fractional exponents (and potentially introduce spurious solutions)
Eq2=16*Eq ; print("Eq2 = ", Eq2)
# Transform log expressions
w0, w1, w2=(SR.wild(u) for u in range(3))
R1={w0*log(w1)+w2:log(w1^w0)+w2}
Eq3=Eq2.subs(R1) ; print("Eq3 = ", Eq3)
# Exponentiate both sides
Eq4=Eq3.operator()(*map(exp, Eq3.operands())).expand_log()
print("Eq4 = ", Eq4)
# Roots = candidate solutions
Cands=Eq4.subs({Cent:100, K:35000/4}).roots(ring=QQbar)
print("Cands = ", Cands)
# Eliminate spurious solutions
Sols=[u[0] for u in Cands
if (C-D).subs({Cent:100, K:35000/4}).subs(y==u[0]).abs()<1e-10]
print("Sols = ", Sols)
# Numerical check
print("Numerical solution : ",
(C-D).subs(Cent=100,K=35000/4).find_root(0,100))


prints :

Eq =  3/16*log(Cent + K) + 5/16*log(Cent) == 1/4*log(Cent + K - y) + 1/4*log(Cent - y)
Eq2 =  3*log(Cent + K) + 5*log(Cent) == 4*log(Cent + K - y) + 4*log(Cent - y)
Eq3 =  log((Cent + K)^3) + 5*log(Cent) == log((Cent + K - y)^4) + 4*log(Cent - y)
Eq4 =  (Cent + K)^3*Cent^5 == (Cent + K - y)^4*(Cent - y)^4
Cands =  [(67.14724076106359?, 1), (133.1013308841840?, 1), (8816.89866911582?, 1), (8882.85275923894?, 1), (99.87573182854424? - 32.97517161162810?*I, 1), (99.87573182854424? + 32.97517161162810?*I, 1), (8850.12426817146? - 32.97517161162810?*I, 1), (8850.12426817146? + 32.97517161162810?*I, 1)]
Sols =  [67.14724076106359?]
Numerical solution :  67.14724076106359


EDIT 4 : In the previous edit, the fragment:

# Get rid of fractional exponents (and potentially introduce spurious solutions)
Eq2=16*Eq ; print("Eq2 = ", Eq2)
# Transform log expressions
w0, w1, w2=(SR.wild(u) for u in range(3))
R1={w0*log(w1)+w2:log(w1^w0)+w2}
Eq3=Eq2.subs(R1) ; print("Eq3 = ", Eq3)


can be replaced by the shorter(but less obvious) :

sage: import sympy
sage: sympy.logcombine((Eq*16)._sympy_())._sage_()
3*log(Cent + K) + 5*log(Cent) == 4*log(Cent + K - y) + 4*log(Cent - y)


This might be a useful addition to Sage's repertoire... (to be considered).

EDIT 5: (ever heard about "comic of repetition" ?...)

An even shorter shortcut :

sage: (16*Eq).maxima_methods().ratsimp()
5*log(Cent) + 3*log(1/4*K + 100) == 4*log(Cent - y) + 4*log(1/4*K - y)


HTH,

more

For some reason, the solution you found is not exactly the same @CyrilleP found, which is also the one found in find_root and sympy's solveset.

( 2021-12-07 22:03:13 +0200 )edit
1

Agreed. See the third and fourth edits.

( 2021-12-08 11:59:55 +0200 )edit