Matlab: Solving a logarithmic equation

Question

I have the following equation that I want to solve with respect to a:

x = (a-b-c+d)/log((a-b)/(c-d))

where x, b, c, and d are known. I used Wolfram Alpha to solve the equation, and the result is:

a = b-x*W(-((c-d)*exp(d/x-c/x))/x)

where W is the is the product log function (Lambert W function). It might be easier to see it at the Wolfram Alpha page.

I used the Matlab's built-in lambertW function to solve the equation. This is rather slow, and is the bottleneck in my script. Is there another, quicker, way to do this? It doesn't have to be accurate down to the 10th decimal place.

EDIT: I had no idea that this equation is so hard to solve. Here is a picture illustrating my problem. The temperatures b-d plus LMTD varies in each time step, but are known. Heat is transferred from red line (CO2) to blue line (water). I need to find temperature "a". I didn't know that this was so hard to calculate! :P

Answer 1

Another option is based on the simpler Wright ω function:

a = b - x.*wrightOmega(log(-(c-d)./x) - (c-d)./x);

provided that d ~= c + x.*wrightOmega(log(-(c-d)./x) - (c-d)./x) (i.e., d ~= c+b-a, x is 0/0 in this case). This is equivalent to the principal branch of the Lambert W function, W₀, which I think is the solution branch you want.

Just as with lambertW, there's a wrightOmega function in the Symbolic Math toolbox. Unfortunately, this will probably also be slow for a large number of inputs. However, you can use my wrightOmegaq on GitHub for complex-valued floating-point (double- or single-precison) inputs. The function is more accurate, fully-vectorized, and can be three to four orders of magnitude faster than using the built-in wrightOmega for floating-point inputs.

For those interested, wrightOmegaq is based on this excellent paper:

Piers W. Lawrence, Robert M. Corless, and David J. Jeffrey, "Algorithm 917: Complex Double-Precision Evaluation of the Wright omega Function," ACM Transactions on Mathematical Software, Vol. 38, No. 3, Article 20, pp. 1-17, Apr. 2012.

This algorithm goes beyond the cubic convergence of the Halley's method used in Cleve Moler's Lambert_W and uses a root-finding method with fourth-order convergence (Fritsch, Shafer, & Crowley, 1973) to converge in no more than two iterations.

Also, to further speed up Moler's Lambert_W using series expansions, see my answer at Math.StackExchange.