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

  • 0
    You could try to solve the equation numerically at a few values, like $y = 0,0.1,0.2,...,35$ and interpolate in between (piecewise or quadratically if it must be continuously differentiable for instance)2011-04-07

2 Answers 2

3

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.

  • 0
    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.2011-04-10
1

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.

  • 1
    @Joel Barsotti: I think a little study of [Galois's biography](http://en.wikipedia.org/wiki/%C3%89variste_Galois) might be in order here.2011-04-07