First time here? Check out the FAQ!

Ask Your Question
0

Solve return an implicit equation

asked 3 years ago

CyrilleP gravatar image

updated 3 years ago

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 ?

Preview: (hide)

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 ( 3 years ago )

1 Answer

Sort by » oldest newest most voted
1

answered 3 years ago

Emmanuel Charpentier gravatar image

updated 3 years ago

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,

Preview: (hide)
link

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 ( 3 years ago )
1

Agreed. See the third and fourth edits.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 3 years ago )

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: 3 years ago

Seen: 356 times

Last updated: Dec 10 '21