Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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

share|improve this question
1  
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
1  
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

2 Answers 2

up vote 3 down vote accepted

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.

share|improve this answer
    
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. –  J. M. 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.

share|improve this answer
    
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. –  J. M. 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
1  
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. –  J. M. Apr 7 '11 at 6:08
1  
@Joel Barsotti: I think a little study of Galois's biography might be in order here. –  Henry Apr 7 '11 at 17:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.