Ask Your Question
0

Solve return an implicit equation

asked 2021-12-06 18:13:58 +0100

CyrilleP gravatar image

updated 2021-12-06 18:14:45 +0100

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 flag offensive close merge delete

Comments

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.

Max Alekseyev gravatar imageMax Alekseyev ( 2021-12-06 18:29:48 +0100 )edit

1 Answer

Sort by ยป oldest newest most voted
1

answered 2021-12-07 17:14:26 +0100

Emmanuel Charpentier gravatar image

updated 2021-12-10 17:32:46 +0100

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,

edit flag offensive delete link more

Comments

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.

cav_rt gravatar imagecav_rt ( 2021-12-07 22:03:13 +0100 )edit
1

Agreed. See the third and fourth edits.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2021-12-08 11:59:55 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2021-12-06 18:13:58 +0100

Seen: 336 times

Last updated: Dec 10 '21