7

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 enter image description here

ROLF
  • 284
  • 2
  • 14
  • 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? – knedlsepp Feb 24 '15 at 17:55
  • 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. – ROLF Feb 24 '15 at 18:13
  • 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. – eigenchris Feb 24 '15 at 18:33
  • 1
    @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. – knedlsepp Feb 24 '15 at 18:48
  • 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: http://octave-specfun.sourcearchive.com/documentation/1.0.9-1/lambertw_8m-source.html Did not benchmark the code, but with a for loop with only 10 iterations and no other loops it should be fast. – Daniel Feb 24 '15 at 19:01
  • 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. – horchler Feb 25 '15 at 02:18
  • @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). – horchler Feb 25 '15 at 16:40
  • 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. – horchler Feb 25 '15 at 16:52

3 Answers3

4

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, W0, 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.

Community
  • 1
  • 1
horchler
  • 18,384
  • 4
  • 37
  • 73
  • 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. – knedlsepp Feb 25 '15 at 01:18
  • 1
    @knedlsepp: Halley's method is a *second order* [Householder class method](http://en.wikipedia.org/wiki/Householder%27s_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)`. – horchler Feb 25 '15 at 02:06
  • 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.) – knedlsepp Feb 25 '15 at 11:13
  • @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`. – horchler Feb 25 '15 at 14:40
  • I still don't get what you mean by *third-order root-finding method* if not the order of convergence? – knedlsepp Feb 25 '15 at 16:47
  • 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. – horchler Feb 25 '15 at 17:06
  • http://s3.amazonaws.com/researchcompendia_prod/articles/e77ac909e9346e07d855c6d006412881-2013-12-23-02-10-26/a20-lawrence.pdf Look for 'order' if you like... Also wikipedia is in line with this: http://en.wikipedia.org/wiki/Root-finding_algorithm – knedlsepp Feb 25 '15 at 17:08
  • To me it just seems you don't like admitting when you are wrong. – knedlsepp Feb 25 '15 at 17:10
  • @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. – horchler Feb 25 '15 at 17:28
  • 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. – horchler Feb 25 '15 at 17:31
  • 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. – knedlsepp Feb 25 '15 at 17:36
  • @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. – horchler Feb 25 '15 at 18:14
  • @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)? – ROLF Mar 09 '15 at 16:31
  • 1
    @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. – horchler Mar 09 '15 at 18:20
  • Thank you! :) I do not need more than one or two digit accuracy. Can your formula above be changed in order to improve accuracy? – ROLF Mar 09 '15 at 20:12
  • @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. – horchler Mar 09 '15 at 21:17
  • @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 – horchler Mar 09 '15 at 21:49
1

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.
knedlsepp
  • 6,065
  • 3
  • 20
  • 41
  • 1
    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. – TroyHaskin Feb 24 '15 at 19:44
  • @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. – knedlsepp Feb 24 '15 at 19:56
  • @ROLF: I don't know what your updated question should tell me? – knedlsepp Feb 25 '15 at 11:16
0

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

GuestF
  • 1