CODEWITHSUNDEEP
matlabsolverexponential
user2349373
user2349373

Reputation: 51

What is the fastest method for solving exp(ax)-ax+c=0 for x in matlab

What is the least computational time consuming way to solve in Matlab the equation:

exp(ax)-ax+c=0

where a and c are constants and x is the value I'm trying to find? Currently I am using the in built solver function, and I know the solution is single valued, but it is just taking longer than I would like.

Upvotes: 3

Views: 1059

Answers (2)

horchler
horchler

Reputation: 18484

Solving your system symbolically

syms a c x;
fx0 = solve(exp(a*x)-a*x+c==0,x)

which results in

fx0 =

   (c - lambertw(0, -exp(c)))/a

As @woodchips pointed out, the Lambert W function has two primary branches, W0 and W−1. The solution given is with respect to the upper (or principal) branch, denoted W0, your equation actually has an infinite number of complex solutions for Wk (the W0 and W−1 solutions are real if c is in [−∞, 0]). In Matlab, lambertw is only implemented for symbolic inputs and thus is very slow method of solving your equation if you're interested in numerical (double precision) solutions.

If you wish to solve such equations numerically in an efficient manner, you might look at Corless, et al. 1996. But, as long as your parameter c is in [−∞, 0], i.e., -exp(c) in [−1/e, 0] and you're interested in the W0 branch, you can use the Matlab code that I wrote to answer a similar question at Math.StackExchange. This code should be much much more efficient that using a naïve approach with fzero.

If your values of c are not in [−∞, 0] or you want the solution corresponding to a different branch, then your solution may be complex-valued and you won't be able to use the simple code I linked to above. In that case, you can more fully implement the function by reading the Corless, et al. 1996 paper or you can try converting the Lambert W to a Wright ω function: W0(z) = ω(log(z)), W−1(z) = ω(log(z)−2πi). In your case, using Matlab's wrightOmega, the W0 branch corresponds to:

fx0 =

   (c - wrightOmega(log(-exp(c))))/a

and the W−1 branch to:

fxm1 =

   (c - wrightOmega(log(-exp(c))-2*sym(pi)*1i))/a

If c is real, then the above reduces to

fx0 =

   (c - wrightOmega(c+sym(pi)*1i))/a

and 

fxm1 =

   (c - wrightOmega(c-sym(pi)*1i))/a

Matlab's wrightOmega function is also symbolic only, but I have written a double precision implementation (based on Lawrence, et al. 2012) that you can find on my GitHub here and that is 3+ orders of magnitude faster than evaluating the function symbolically. As your problem is technically in terms of a Lambert W, it may be more efficient, and possibly more numerically accurate, to implement that more complicated function for the regime of interest (this is due to the log transformation and the extra evaluation of a complex log). But feel free to test.

Upvotes: 1

user85109
user85109

Reputation:

Just wanting something to run more quickly is insufficient for that to happen.

And, sorry, but if fzero is not fast enough then you won't do much better for a general root finding tool.

If you aren't using fzero, then why not? After all, that IS the built-in solver you did not name. (BE EXPLICIT! Otherwise we must guess.) Perhaps you are using solve, from the symbolic toolbox. It will be more slow, since it is a symbolic tool.

Having said the above, I might point out that you might be able to improve by recognizing that this is really a problem with a single parameter, c. That is, transform the problem to solving

exp(y) - y + c = 0

where

y = ax

Once you know the value of y, divide by a to get x.

Of course, this way of looking at the problem makes it obvious that you have made an incorrect statement, that the solution is single valued. There are TWO solutions for any negative value of c less than -1. When c = -1, the solution is unique, and for c greater than -1, no solutions exist in real numbers. (If you allow complex results, then there will be solutions there too.)

So if you MUST solve the above problem frequently and fzero was inadequate, then I would consider a spline model, where I had precomputed solutions to the problem for a sufficient number of distinct values of c. Interpolate that spline model to get a predicted value of y for any c.

If I needed more accuracy, I might take a single Newton step from that point.

In the event that you can use the Lambert W function, then solve actually does give us a solution for the general problem. (As you see, I am just guessing what you are trying to solve this with, and what are your goals. Explicit questions help the person trying to help you.)

solve('exp(y) - y + c')
ans =
c - lambertw(0, -exp(c))

The zero first argument to lambertw yields the negative solution. In fact, we can use lambertw to give us both the positive and negative real solutions for any c no larger than -1.

X = @(c) c - lambertw([0 -1],-exp(c));
X(-1.1)
ans =
     -0.48318      0.41622

X(-2)
ans =
      -1.8414       1.1462

Upvotes: 3

Related Questions