This question on StackOverflow has led me to ask myself what would be involved in solving an equation like this one using tools available to the Python programmer. $$ \frac{0.125567841}{d^{2.25}} = \frac{2.513274d+0.10053+2.11\pi}{d+0.04}$$

What are the fundamental numerical tools that one would use to get the one real and four complex roots of this equation (WolframAlpha)? I have looked (hard) at the Mathematica documentation but cannot discern what they might be. It's probably mainly ignorance on my part.

  • $\begingroup$ Numerically or symbolically? $\endgroup$
    – rubenvb
    Oct 23 '17 at 22:40
  • $\begingroup$ @rubenvb: Either way. Thanks for responding. $\endgroup$
    – Bill Bell
    Oct 23 '17 at 22:42
  • $\begingroup$ Does Mathematica solve this symbolically? $\endgroup$
    – nicoguaro
    Oct 23 '17 at 23:23
  • $\begingroup$ Mathematica is a proprietary closed-source software package, so there is no documentation that is freely available about "how" Mathematica does most things. There is abundant literature on root-finding methods for polynomials, both symbolically and numerically, and study of such algorithms will likely prove to be a rewarding project. However as a "black box" we cannot be sure at a low level how Mathematica accomplishes either of these. $\endgroup$
    – hardmath
    Oct 24 '17 at 1:49
  • $\begingroup$ @nicoguaro: In case you're asking me, I don't know. $\endgroup$
    – Bill Bell
    Oct 24 '17 at 3:29

You can rewrite with integer exponents through a change of variable. This looks like $d=x^r$ where $r$ is the LCM of the denominator of the rationalized exponents. In your case we have $r=4$ and:

$$6.72427 - 0.125568 x^4 + 2.51327 x^{13}=0$$

Which can be solved in the usual ways. In Mathematica there is a clear distinction between symbolic and numerical solutions. You have floating point numbers bouncing around, so the solution that pops out is a result of numerical approximation. These are described in their extended internal documentation: http://reference.wolfram.com/language/tutorial/SomeNotesOnInternalImplementation.html#2118

The 13 roots lie along a circle in the complex plane: enter image description here


As the other Answer already touches on the possibility of a symbolic root-solver being applied to this particular equation (by transforming into a polynomial form, albeit of degree $\ge 5$), I'll make some remarks about numerical root-finding.

When we want to find real roots of a continuous function, the Intermediate Value Theorem is a basic tool. On a closed bounded interval where the real function is continuous and changes sign between endpoints, there will be a real root somewhere in the interior of that interval. The Bisection Method allows us to approximate at least one such root in that interval.

Here we can transform the equation into finding roots of a real continuous function by cross-multiplying and collecting terms on one side:

$$ \frac{0.125567841}{d^{2.25}} = \frac{2.513274d+0.10053+2.11\pi}{d+0.04} $$

$$ 0.125567841(d+0.04) = (2.513274d+0.10053+2.11\pi)d^{2.25} $$

$$ F(d) \equiv 2.513274d^{3.25} + (0.10053+2.11\pi)d^{2.25} - 0.125567841d - 0.00502271364 = 0 $$

This form of the problem is easier to analyze, but we should check whether the cross-multiplication introduced artifact roots, i.e. roots for the revised equation that do not satisfy the original equation. When we multiply both sides by the denominators $d^{2.25}$ and $d+0.04$, it is possible that this introduces one or more new roots where such factors are zero, i.e. where $d=0$ or $d=-0.04$. Fortunately $F(0)\neq 0$ and $F(-0.04)\neq 0$, so we have not introduced artifact roots in this case.

One should be careful to define what (if anything) $d^{2.25}$ means for negative real values of $d$. Inevitably one must abandon any real value for this expression because it requires taking square roots of a negative number. This is a particular weakness of relying on tools such as Wolfram Alpha, because it attempts to be helpful by aggressively assuming that some sense is to be made of the expression even for negative arguments $d$.

Fortunately in this case what roots Wolfram Alpha finds are not on the negative half of the real axis for $d$. There is a positive real root $d\approx 0.06087$, and the other roots found are in complex conjugate pairs.

The real root can be found by standard methods such as bisection once we realize that $F(d)$ is continuous and changes sign on the interval $[0,1]$. That is, $F(0)=-0.00502271364$ and $F(1)\approx 9.11197$, so there must be an intermediate value of $d_0$ where the function $F(d_0)=0$. Once the location of the root $d_0$ is suitably narrowed by bisection, etc., more rapidly converging methods (like Newton's method) can be used.

When a function like $F(d)$ with real coefficients on rational powers has complex roots, these will occur naturally in complex conjugate pairs. One can then replace the one dimensional search for a real root with a two dimensional search for a complex conjugate root pair $d = x \pm iy$. A classic approach to this is Bairstow's Method, which uses Newton-Raphson iterations to locate a complex conjugate root pair using only real arithmetic.

Just as Python users interested in symbolic root-finding should be familiar with sympy, those interested in numerical root-finding will be attracted to the numpy library. In many problems a synthesis of the two approaches is advantageous, e.g. using sympy to get derivatives of a function which can be evaluated by routines in numpy.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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