0

I have an equation,

(1 - y)/y^4 = x^2

and I need to plot y(x). I have solved it using this snippet:

from sympy import symbols, Eq, solve

x, y = symbols('x y')
equation = Eq((y-1)*y**(-4), x**2)

solution = solve(equation, y)
print(solution)

And I get this monolith of an output:

[Piecewise((-sqrt(-(-1/x**4)**(1/3))/2 - sqrt((-1/x**4)**(1/3) - 2/(x**2*sqrt(-(-1/x**4)**(1/3))))/2, Eq(x**(-2), 0)), (-sqrt(2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))/2 - sqrt(-2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) - 2/(x**2*sqrt(2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))) - 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))/2, True)), Piecewise((-sqrt(-(-1/x**4)**(1/3))/2 + sqrt((-1/x**4)**(1/3) - 2/(x**2*sqrt(-(-1/x**4)**(1/3))))/2, Eq(x**(-2), 0)), (-sqrt(2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))/2 + sqrt(-2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) - 2/(x**2*sqrt(2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))) - 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))/2, True)), Piecewise((sqrt(-(-1/x**4)**(1/3))/2 - sqrt((-1/x**4)**(1/3) + 2/(x**2*sqrt(-(-1/x**4)**(1/3))))/2, Eq(x**(-2), 0)), (sqrt(2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))/2 - sqrt(-2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(x**2*sqrt(2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))) - 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))/2, True)), Piecewise((sqrt(-(-1/x**4)**(1/3))/2 + sqrt((-1/x**4)**(1/3) + 2/(x**2*sqrt(-(-1/x**4)**(1/3))))/2, Eq(x**(-2), 0)), (sqrt(2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))/2 + sqrt(-2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(x**2*sqrt(2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3) + 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))) - 2/(3*x**2*(sqrt(-1/(27*x**6) + 1/(256*x**8)) + 1/(16*x**4))**(1/3)))/2, True))]

I would like to plot it over the bounds 0 <= y <= 1, and 0 <= x <= some_large_number. I haven't been able to figure out how to do it, but I'm sure it's got to be easy for me, if a PITA for my CPU.

Also, it is of no interest to me to solve this analytically. I would be only too delighted to do it numerically, but I don't know how to do that either. Sympy was just a solution that presented itself to me.

Frank Harris
  • 587
  • 2
  • 16

1 Answers1

0

Sympy seems to be mixing in complex solutions. I'm unsure how to cope with them in this case. Note that for y smaller than 1, the LHS is negative, so never equal to the square of a real number.

If, however, you are mostly interested in a plot, plot_implicit can come in handy:

from sympy import symbols, Eq
from sympy.plotting import plot_implicit

x, y = symbols('x y') # symbols('x y', real=True, positive=True) doesn't seem to make a difference
equation = Eq(((y - 1) * y ** (-4)), x ** 2)

plot_implicit(equation)

example plot

PS: To find the rightmost point of the curve, differentiating the left-hand side solve(equation.lhs.diff()) gives the solution y = 4/3, corresponding to x = 2*sqrt(3)/3. Finding numerical solutions for y >= 1 is straightforward. nsolve can find a numerical solution for a given x0 smaller than 2*sqrt(3)/3, e.g. nsolve(equation.subs(x, .1), 1) gives 1.01042350456793.

JohanC
  • 71,591
  • 8
  • 33
  • 66