The goal is to find all real valued solutions to the form of any equation of the form
where all and all (Next step would be to extend it to handle complex numbers or complex conjugates.) Note that is one of the B terms is 0, then the term is constant.These equations frequently arise in physics from solutions of systems of differential equations. In my case I needed to solve the to simulate a heat exchanger.
As far as I know, the method used has not been described elsewhere, so I will christen it "Eng's method" after yours truly.
The basic method is as follows:
- Sort the terms by ascending value of exponent:
- Find a range where there could possibly be a root. To do so, consider that as x increases, the Nth term will begin growing faster than all other terms, so find a value of x where Actually we will do the following. Count the number of terms with an opposite sign as . Call this number P. For all P terms calculate So basically at we have guaranteed that the fastest growing term is growing P times faster than any other term with opposite sign. So there will definitely not be a root for .
- Using the same sort of reasoning, looking at the slowest growing term and find an x where this term is Q times more than all other terms with opposite sign. The minimum such term is . Since the slowest growing term is also the slowest to shrink as we move towards , we can confidentially say that for the sign of the first term will dominate and we will never cross .
- Take repeated derivatives (k of them) of .
-
At this point we know that for the kth derivative a single term is dominant over the entire range where there could possibly be a root. What this means is that the sign of the kth derivative is either always positive or always negative for . Therefore it follows that the (k-1)th derivative must be monatomic over this range. This is great news because for a monatomic function we can find numerically find roots using algorithms like bisection of Brent's method. In fact, all we have to do is evaluate the (k-1)th derivative at and at . If these endpoints are the same sign, then we know that the (k-1)th derivative is also the same sign over the entire range. However, if they are the same sign, then we can use bisection to find the root of the (k-1)th derivative with guaranteed convergence to any desired degree of accuracy. If we can find where then this will leave us with two range, and .
-
Now we either know there is one range where the (k-1)th derivative has a consistent sign or two ranges where the (k-1)th derivative has a consistent sign for the whole range (but one range is positive and the other is negative). Either way, we know that on these one or two ranges we can use bisection (or another bracketed root finding method) to find if and where .
-
We continue in this manner, breaking our original range into smaller ranges when we find that the current Kth derivative is 0 and then proceeding to use bracketed root finding at one less level derivative. At some point, we will keep backing out of layers of derivatives until we are just doing bracketed root finding on the original function.
TLDR;
- Find a range where roots could possibly occur.
- Take a bunch (N) derivatives until its obvious that one term is dominant in the range where roots occur.
- With exponents this always happens eventually.
- This tells us that the (N-1)th derivative is a monatomic function so we can find a root (if it has one) via bisection.
- Now we have some more regions on which we know the the (N-2) derivative is a monatomic function.
- Apply logic recursively until we have regions of the original curve where we can find roots using bisection. This way we won't miss any roots.
Sum of exponent functions are defined as arrays of terms:
[{A: 1, B: 1}]
would correspond to and [{A: 1, B: 1},{A: -2, B: 0}]
corresponds to
Two functions are provided, evaluateExponentFormula
and findRootOfExponents
evaluateExponentFormula(formula, x)
takes a formula returns it evaluated it at x.
findRootOfExponents(xMin, xMax, formula)
takes a formula and returns all the roots between xMin and xMax (you could use -Infinity and Infinity to get them all).
> node solve-exponentials-tests.mjs
Warning: at this point, this is a proof of concept and provided as is. I am sure there are many unhandled edge cases.