3
$\begingroup$

$$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

$\endgroup$
8
  • 1
    $\begingroup$ 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. $\endgroup$ Apr 6, 2011 at 23:17
  • $\begingroup$ @Joel: do you mean to solve for $x$? Use numerical methods. $\endgroup$ Apr 6, 2011 at 23:17
  • $\begingroup$ Yes I need to solve for x. $\endgroup$ Apr 6, 2011 at 23:46
  • 1
    $\begingroup$ bit.ly/hRsmVM (just a direct link to a Wolfram Alpha page with the answer.) $\endgroup$
    – Alon Amit
    Apr 6, 2011 at 23:52
  • $\begingroup$ 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. $\endgroup$ Apr 7, 2011 at 0:07

2 Answers 2

3
$\begingroup$

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"}]

Residual plot

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.

$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – Myself
    Apr 8, 2011 at 18:01
  • $\begingroup$ @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. $\endgroup$
    – whuber
    Apr 8, 2011 at 19:15
  • $\begingroup$ @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. $\endgroup$
    – whuber
    Apr 8, 2011 at 19:18
  • $\begingroup$ 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. $\endgroup$ Apr 10, 2011 at 18:00
1
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$ Apr 7, 2011 at 2:41
  • $\begingroup$ What I need to be able to do is find a given X for a Y value. $\endgroup$ Apr 7, 2011 at 2:45
  • 1
    $\begingroup$ 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 :) $\endgroup$
    – GWu
    Apr 7, 2011 at 5:12
  • $\begingroup$ 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. $\endgroup$ Apr 7, 2011 at 6:08
  • 1
    $\begingroup$ @Joel Barsotti: I think a little study of Galois's biography might be in order here. $\endgroup$
    – Henry
    Apr 7, 2011 at 17:26

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .