Solve return an implicit equation

edit

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,