Solving an algebraic equation with sympy

Question

I have this algebraic equation that can be solved in WolframAlpha:

but it cannot be solved in Python.

There is this error:

NotImplementedError: multiple generators [10**x, 2**x]
No algorithms are implemented to solve equation 
       (1/2)**(1 - x) - 10**(x - 1)/4 - 3/4

I was wondering what is wrong with my code?

x = symbols('x')
eq1 = Eq(0.5**(1-x)-1-((1/4)*(0.1**(1-x)-1)),0)
sol = solve((eq1),(x))
sol

Answer 1

When solve fails, you can attempt to find a numerical solution with nsolve. First, you need to find reasonable initial guesses. In this case a simple plot is enough:

from sympy import *
x = symbols("x")
eq = 0.5**(1-x)-1-((1/4)*(0.1**(1-x)-1))
plot(eq, ylim=(-5, 5))

It appears there is at least one root at x=1. Let's zoom in:

plot(eq, (x, 0, 2))

There are two roots. Let's find them:

initial_guesses = [1, 1.5]
r = [nsolve(eq, x, i) for i in initial_guesses]
print(r)
# out: [1.00000000000000, 1.21889091552120]

Answer 2

This is not really an algebraic equation... And the solutions of equations like that are not guaranteed to be algebraic numbers. Said differently, apart from the evident (and exact) solution x = 1 any other solution is likely to be a transcendental number.

As this equation has no nice property (at least I could not find any at first sight), I think that it cannot be solved analytically, and SymPy is correct at saying that it cannot solve it: it has no algorithm for that... because no is known by mathematicians!

I do not know WolframAlpha, but I would assume that it just solved it numerically. And it is possible in Python too, either by hand with a little dichotomy followed by the Euler method, or with sympy.nsolve as shown in @Davide_sd's answer, or with scipy, or...

Answer 3

You can write your own solver to solve the equation by single-point iteration.

It's helpful to rewrite the equation first:

0.5 * 2^x -3/4 -0.025 * 10^x = 0

or, multiplying by 40 and rearranging,

20 * 2^x - 30 = 10^x

You can rearrange to make either x-containing term the subject and then take logs (any base) to give two possible iterative formulae:

x = log(20 * 2^x - 30)/log(10)

or

x = log(1.5+0.05 * 10^x)/log(2)

The first is appropriate for single-point iteration when x>=1, the latter when x<=1. Moreover, since 10^x is positive then 20 * 2^x - 30 must be also, which means that x > log(1.5)/log(2).

In the following code enter a starting point less than 1 to get the smaller solution and a starting point greater than 1 to get the second.

from math import log

def RHS( x ):
    if x > 1.0:
       return log( 20 * 2 ** x - 30 ) / log( 10.0 )
    else:
       return log( 1.5 + 0.05 * 10 ** x ) / log( 2.0 )

TOLERANCE = 1.0e-20;
XMIN = log( 1.5 ) / log( 2.0 )

print( "Enter a starting point for iteration (greater than ", XMIN, "): " )
x = float( input() )
xprev = x + 1.0

while abs( x - xprev ) > TOLERANCE:
    xprev = x
    x = RHS( x )

print( "Solution: {:13.6f}".format( x ) )

Output examples:

Enter a starting point for iteration (greater than  0.5849625007211562 ): 
0.7
Solution:      1.000000


Enter a starting point for iteration (greater than  0.5849625007211562 ): 
10
Solution:      1.218891