# Help with a quintic polynomial

$$y = 0.10 + 4.060264x - 6.226862x^2 + 48.145864x^3 - 60.928632x^4 + 49.848766x^5$$

I need to be able to solve this equation for $x$.

I've looked around and seem to be failing miserably and solving this myself.

I'll have a $y$ value (likey between 0-35) and I need to find an exact $x$ (likely between 0-1).

Thanks

-
Do you mean you need to solve y=0 for an x between 0 and 1? Wolfram Alpha says there is only one real root, and it is negative. –  Matthew Conroy Apr 6 '11 at 23:17
@Joel: do you mean to solve for $x$? Use numerical methods. –  Weltschmerz Apr 6 '11 at 23:17
Yes I need to solve for x. –  Joel Barsotti Apr 6 '11 at 23:46
bit.ly/hRsmVM (just a direct link to a Wolfram Alpha page with the answer.) –  Alon Amit Apr 6 '11 at 23:52
sorry I guess I'm not being able to clearly describe my problem. I need to be able to solve it given a Y value between 0 and 35. –  Joel Barsotti Apr 7 '11 at 0:07

One way to obtain a "symbolic version," as requested in the comments, is to compute some relatively simple approximation and polish it with a Newton-Raphson step. Because this function is smooth and monotonic for $0 \le y \le 35$ this is going to work very well.

In fact, a least-squares fit of the functional form $a \log(b + c(y+1)^{1/5} + d(y-e)^2$ to the solutions for $y=0, 1, \ldots, 35$ already gets close: most of the errors are less than 0.0003 . One Newton-Raphson step is a rational function of this expression of degree 5 (numerator) and 4 (denominator), thereby expressible in terms of 11 parameters derived from the original polynomial. The residuals of this 16-parameter expression range from $-6 10^{-6}$ to $2.7 10^{-7}$, which is close to the precision of the original polynomial coefficients. For $y \ge 4$ the errors are all less than $10^{-7}$, which is as good as one can hope for.

To find this solution in Mathematica, begin by generating the array of solutions for $y=0, 1, \ldots, 35$:

Clear[x, y];
roots = x /. Table[FindRoot[-y + 0.10 + 4.060264 x - 6.226862 x^2 +
48.145864 x^3 - 60.928632 x^4 + 49.848766 x^5, {x, .5}], {y, 0, 35}]


Fit the initial simple model (using some eyeball guesses for the parameters):

Clear[a, b, c];
model = a Log[b + c  y^(1/5)] +  d (y - e)^2;
fit = FindFit[roots, model, {{a, .5}, {b, 1}, {c, .1}, {e, 18}, {d, .0001}}, y]


Create a Newton-Raphson step for a function f at the argument a:

Clear[nr];
nr[f_, a_] := (x - f[x]/D[f[x], x]) /. x -> a


Use it to improve the model:

Clear[x];
x[z_] := ( nr[f[#] - y + 1 &, model /. fit ]) /. y -> (z + 1)


(The shift to y-1 from y is needed because Mathematica starts indexing at 1, not 0.) The model works well for $1 \le y \le 35$ and exceptionally well for $y \ge 4$.

g = Table[x[y], {y, 1, 36}];
ListPlot[roots - g, PlotRange -> {Full, Full},
PlotStyle -> PointSize[0.015], DataRange -> {0, 35},
AxesLabel -> {"y", "Error"}]


If you need better solutions for $y \lt 4$, you could similarly fit a simple model plus a Newton-Raphson polish to this range of values alone.

-
This is a nice way to deal with this problem, I think! But I'm wondering: where does the form $a \log(b + c(y+1)^{1/5} + d(y-e)^2$ come from? –  Myself Apr 8 '11 at 18:01
@My It can be anything that works. When you plot the inverse of the polynomial it looks roughly logarithmic for $0 \lt x \lt 1$. Including the 1/5 power is a naive attempt to deal with the degree of the polynomial--but it works. The +1 is there to avoid problems with negative arguments. After fitting an initial model with $a, b, c,$ the residuals for the range $4 \lt y \lt 35$ describe almost a perfect parabola, which is fit with $d$ and $e$. That gave 3-4 digits of accuracy. Because each NR step usually doubles the sig figs, 3-4 is just good enough for one polishing step to succeed. –  whuber Apr 8 '11 at 19:15
@My You also found a slight inconsistency in the text: I gave up on the $(y+1)^{1/5}$ part and reverted to $y^{1/5}$ (shown in the code blocks) because the latter fit better and was less complicated. –  whuber Apr 8 '11 at 19:18
Nice! I did something similar, but I constructed a Padé approximant instead of doing least squares. I wound up with $$\frac{\frac1{562}(y-10)^2+\frac3{35}(y-10)+\frac{16}{23}}{\frac1{840}(y-10)^2+‌​\frac3{34}(y-10)+1}$$ ... a bit rough around the ends, but quite good in the middle. –  Ｊ. Ｍ. Apr 10 '11 at 18:00

It is a quintic, and so must have at least one real root. The coefficient of $x^5$ is positive so $y$ is negative if $x$ is large and negative and $y$ is positive if $x$ is large and positive.

"By inspection" you will get negative $y$ for $x=-2$ and you obviously get a positive $y$ for $x=0$. In fact it passes through the points $(-2, -2988.113512)$ and $(0, 0.1)$, so there is a root between them so try halfway, finding the point $(-1,-169.110388)$ and continue halving the interval, so next $(-0.5,-14.8708939375)$; you don't need such precision, and in fact only need the sign of $y$. Just keep halving between two values of $x$ which give opposite signs for $y$ until you get an answer which is precise enough for you.

There are faster root-finding algorithms, but this bisection method will work whenever you have a continuous function and two points with $y$ having opposite signs.

-
To add: it is profitable to exploit Descartes' rule of signs (alternatively, Fourier-Budan) for localizing the root before entering some iterative scheme. If bisection is too slow, there is a way to properly merge the speed of Newton-Raphson with the sureness of bisection, if one is so inclined. However, if you manage to find a good approximation in the first place, it's probably safe to use Newton-Raphson unfettered. –  Ｊ. Ｍ. Apr 7 '11 at 2:41
What I need to be able to do is find a given X for a Y value. –  Joel Barsotti Apr 7 '11 at 2:45
Then move the given y value to the other side and repeat the same procedure. I hope you're not looking for a closed form formula for x in terms of y, when you said "exact X" in the original post :) –  GWu Apr 7 '11 at 5:12
Another way to say GWu's comment: You do not want the closed form solution here (most certainly if your work is numeric instead of symbolic). It ain't pretty. –  Ｊ. Ｍ. Apr 7 '11 at 6:08
@Joel Barsotti: I think a little study of Galois's biography might be in order here. –  Henry Apr 7 '11 at 17:26