Harry
Harry

Reputation: 99

Finding x with a Python-script in a difficult equation where x is only variable

I have some difficult equations where x is the only variable. For example this one:

atan(12/x)-atan(SQRT(7x)/31)=atan((3x+942+sqrt(2x^2+2))/(4+sqrt(5+x^2)))

In software like WolframAlpha I can paste it in at it will find the x for me easily.

X found to be around 0.00523 by Wolfram Alpha

The only idea I have now is to do some sort of brute force on each side with small iterations until Left side = Right side. I have no idea how to elegantly solve this equation using python.

Does anyone where I should start?

(Disclaimer: I have written some basic python scripts and studied some math, but my math skills and python skills are probably in the bottom 1% in the Stack Overflow-community.)

Best regards,

Upvotes: 1

Views: 673

Answers (1)

coproc
coproc

Reputation: 6257

This is a difficult equation, both symbolically and numerically. Mathematica has very advanced algorithms. nsolve from sympy does not find the root without help, but you can do the following:

from sympy import *
x = symbols('x')
e0 = atan(12/x)-atan(sqrt(7*x)/31)-atan((3*x+942+sqrt(2*x**2+2))/(4+sqrt(5+x**2)))
nsolve(e0,x,0)
# gives ZeroDivisionError
nsolve(e0,x,0.001)
# gives 0.00523293635181032 - 1.6278511388129e-27*I
nsolve(e0,x,0.005)
# gives 0.00523293635181032

So you have to guess the range of the zero rather closely.

To make the search for the root easier, you can get rid of atan by taking the tan of both sides and use tan(a-b) = (tan(a) - tan(b)) / (1 + tan(a)tan(b)) :

t1 = 12 / x
t2 = sqrt(7*x) / 31
t3 = (3*x+942+sqrt(2*x**2+2)) / (4+sqrt(5+x**2))
# atan(t1) - atan(t2) = atan(t3) is equivalent to e1 = 0:
# (note: tan(atan(t1) - atan(t2)) = (t1-t2)/(1+t1*t2) )
e1 = simplify((t1-t2)/(1+t1*t2) - t3)
# gives (sqrt(x)*(31*sqrt(x) + 12*sqrt(7))*(-3*x - sqrt(2*x**2 + 2) - 942) + ...)/(sqrt(x)*...)

Now we have to give nsolve one more hint to make the root search easier, namely to take only the numerator of the expression:

e1_num = e1.as_numer_denom()[0]
# gives (sqrt(x)*(31*sqrt(x) + 12*sqrt(7))*(-3*x - sqrt(2*x**2 + 2) - 942) + (-sqrt(7)*x**(3/2) + 372)*(sqrt(x**2 + 5) + 4)
nsolve(e1_num,x,0)
# gives 0.00523293635181032

Upvotes: 1

Related Questions