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

Comenius

How often and in what context do you need to solve this equation for a? As part of some other solver? Are the values of x,b,c,... which you want to evaluate all already known? – Elianore

I know nothing about the lambert W function. However, I didn't find your answer very helpful, IrrationalPerson. @knedlsepp: I'm using in for solving a heat exchanger in a heat pump. I'm kinda using it in some other solver. I guess the function is used at least 30.000 times, most likely more (I'm doing hourly calculations for a year). As mentioned, x, b, c and d are known. – Comenius

Could you edit your post to provide some realistic values for x,b,c,d? My first naive thought is to try to use fzero() to find the root(s) of (a-b-c+d)/ln((a-b)/(c-d)) - x, but this can lead to complex answers for some of the random parameter values I've tried. – Devland

@ROLF: Have you already tried evaluating lambertw for a vector instead of a single value? Instead of for i = 1:100, lambertw(i), end, you can use lambertw(1:100), which is about 30 times faster than the original. For 30.000 values, this takes about 30 seconds. – Elianore

To add another Idea: lambertw in matlab uses symbolic math, that's a huge overhead. Switch to a numeric implementation. The octave version might be the easiest to port: octave-specfun.sourcearchive.com/documentation/1.0.9-1/… Did not benchmark the code, but with a for loop with only 10 iterations and no other loops it should be fast. – Foramen

Are there any constraints/relationships between a, b, c, d, and/or x? Are there any assumptions on these variables that limit their domain, i.e, do you know if any/all are real-valued or not, are any/all non-negative or or greater/lesser than some value? In some cases, such constraints and assumptions may allow you to simplify things. – Sirenasirenic

@ROLF: Okay, from what I can see, all of the parameters, including x, are real-valued (they're physical values) and positive. Also, a>b, c>d always given the logarithmic mean temperature difference problem. Also, using a function specific to the Lambert W or Wright Omega may be faster and more reliable than using a generic root solver like fzero, but that could work (just be sure to validate it with something else). – Sirenasirenic

Other possibilities to speed things up: 1) try using single precision (this isn't always faster and you have to be careful when and how you convert back and forth to double so as not to slow things down unnecessarily or introduce more imprecision), 2) I don't know if or which of your parameter are scalars or vectors, but it would simplify your equation slightly if you solved x=(delta_ab-delat_cd)/log(delat_ab/delta_cd) for delta_ab and then solved for a from delta_ab=a-b. If this would sped things up or not would depend on your implementation, and the dimensions of your parameters. – Sirenasirenic

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.

Sirenasirenic

I did some testing on the functions to get a better overview: In the interval [-1/e,0] your start values could speed up the process, as Moler's implementation takes about 6 iterations (as opposed to two). In practice though it is already quite fast, as it is fully vectorized. One thing about your answer: I don't know if the algorithm of FSC you mentioned is of even higher order, but Halley's method is a third-order method. – Elianore

@knedlsepp: Halley's method is a second order Householder class method. The rate of convergence is one higher than the order, hence Halley's method has cubic convergence. The Lambert_W function will require more iterations to reach a given tolerance for certain critical values. The example code in my Math.StackExchange answer is just an illustration of using series expansions, not code to be used for comparing speed with another function. The plot shows the more naïve approach uses more iterations, particularly around -1/exp(1). – Sirenasirenic

Oh. the second order Householder class method: Halley's method? This might be confusing, as usually a nth-order method is about the order of convergence. (Also the FSC paper isn't really directly using the Householder method applied to this problem, but provides two algorithms of third and fourth order (of convergence) resp.) – Elianore

@knedlsepp: I agree the order the thing is confusing, but that's not how I learned it. If I want to specify the convergence then I alway say explicitly rate of convergenge or order of convergence to be clear. And, as my answer states, it's the Lawrence, et al. paper that cites the FSC paper as the basis for their iterative algorithm that I implemented in wrightOmegaq. – Sirenasirenic

I still don't get what you mean by third-order root-finding method if not the order of convergence? – Elianore

I'm sorry you don't like it. I provided a link above that says much more than I can fit in a comment here on this: the order relates to the minimum order of continuous differentiability required of the function that one is interested in finding the root of, i.e., the order is intrinsic to the equation that defines the method. – Sirenasirenic

s3.amazonaws.com/researchcompendia_prod/articles/… Look for 'order' if you like... Also wikipedia is in line with this: en.wikipedia.org/wiki/Root-finding_algorithm – Elianore

To me it just seems you don't like admitting when you are wrong. – Elianore

@knedlsepp: I'm trying to be helpful here and clarify my question to someone who is asking questions that aren't even germane to the OP's original one so I don't know how you get off being accusatory. I provided sources and clear explanations. Can you you imagine that you might be wrong or that there is no "right" answer or that you just don't understand what I've written? What is your goal? I think my definition is clear and commonly used. Yes, many may use nebulous language. I really don't know what your problem is or why you keep going on and on about nothing. – Sirenasirenic

The Wikipedia article on root-finding specifically refers to order of convergence, not just order for Newton-like methods, which is exactly what I've said above. This entire discussion is rather pointless. – Sirenasirenic

I think your "definition" of order is nonstandard. Even the paper you mentioned (which I have linked above) refers to the FSC-method as a fourth-order method. – Elianore

@knedlsepp: I don't think it's non standard at all. The order of differentiability is fundamental. And of course the rate/order of convergence is only an upper bound if/when the system actually converges. People can figure out what the writer means in most cases so it's usually not a problem. I'm in academia – all sorts of less than clear stuff gets into papers that should have been caught by reviewers but wasn't. I've updated the wording in my answer to hopefully make everything super explicit. – Sirenasirenic

@horchler: I'm not that into math, so I do not understand the discussion or the math behind your answer. Do I get it right if both wrightOmega and wrightOmegaq can replace the lambertw? In my code, can I simply replace lambertw with one of the two other functions (I understand that wrightOmegaq is not build-in)? – Comenius

@ROLF: You can ignore the comments above -they have nothing to do with your actual question, only the wording in my answer. Best I can tell from your question, the Wright omega function should work fine. I tried some numeric values and modified my formula at the top a bit. You should plug in some of your own values and compare my formula above with the one in your question using lambertw. The two should return very similar results. Matlab's wrightOmega may return a tiny imaginary part due to numeric error (use real to get rid of it). wrightOmegaq doesn't have this issue. – Sirenasirenic

Thank you! :) I do not need more than one or two digit accuracy. Can your formula above be changed in order to improve accuracy? – Comenius

@ROLF: I'm not sure what you mean. I'm guessing you're using double precision floating-point values in Matlab. The output of all or of the methods should be within about machine epsilon (eps) of the true value. If you want less accuracy, then you'll need to use vpa and digits with the symbolic functions lambertw and wrightOmega. For my numeric wrightOmegaq you can try adjusting the value of tol in the code to something larger, like 1e-8 or 1e-4. You may be able to do something similar with Cleve Moler's Lambert_W function too. I don't know how much these will affect speed. – Sirenasirenic

@ROLF: Also, if you provide more exact ranges and relationships for c, d, and x, it may be possible to provide a series expansion solution for this. I assume c>d. Do you know if c-d>0 or c-d>a or c-d>b? What is the range of x? Is x>0 always? Do you know if x>c-d or if x<c-d always? If you look at my code for wrightOmegaq, you'll see that various series expansions are used for initial guesses. These would be good to a few digits of accuracy in most cases, but they depend on the value of the input argument, which can be complex-valued: log(-(c-d)./x)-(c-d)./x. – Sirenasirenic

Answer 2

Two (combinable) options:

• Is your script already vectorized? Evaluate the function for more than a single argument. Executing for i = 1:100, a(i)=lambertw(rhs(i)); end is slower than a=lambertw(rhs).
• If you are dealing with the real valued branch of LambertW (i.e. your arguments are in the interval [-1/e, inf) ), you can use the implementation of Lambert_W submitted by Cleve Moler on the File Exchange.

Elianore

Good catch on the Lambert_W function from Cleve Moler. I was going to recommend a combination of using exp and fzero, but Moler is doing essentially the same thing but with a cubic method, so it's more than likely faster. – Chalfant

@TroyHaskin: Yup, had a quick look at his blog entry on mathworks. It's quite unusual for me to see methods of higher order than newton's iteration being used though. – Elianore

@ROLF: I don't know what your updated question should tell me? – Elianore

Answer 3

Do you know the mass flow rates at both sides of the heat exchanger at each time-step? If yes, temperature 'a' can be solved by the 'effectiveness-NTU' approach which does not need any iteration, rather than the LMTD approach. Reference: e.g. http://ceng.tu.edu.iq/ched/images/lectures/chem-lec/st3/c2/Lec23.pdf

Syzygy