Matlab：求解对数方程答案

【问题标题】：Matlab: Solving a logarithmic equationMatlab：求解对数方程
【发布时间】：2015-02-24 17:48:41
【问题描述】：

关于a，我有以下等式要求解：

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

其中x、b、c 和d 是已知的。我用 Wolfram Alpha 解方程，结果是：

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

其中W 是产品日志函数（Lambert W 函数）。在Wolfram Alpha page 上可能更容易看到它。

我使用 Matlab 内置的lambertW 函数来求解方程。这相当慢，是我脚本中的瓶颈。还有另一种更快的方法吗？它不必精确到小数点后 10 位。

编辑：我不知道这个方程很难解。这是一张说明我的问题的图片。温度 b-d 加上 LMTD 在每个时间步长中都不同，但是是已知的。热量从红线（CO2）转移到蓝线（水）。我需要找到温度“a”。没想到这么难计算！ :P

【问题讨论】：

您需要多久以及在什么情况下为a 求解这个方程？作为其他求解器的一部分？您要评估的x,b,c,... 的值是否都已知？
我对 Lambert W 函数一无所知。但是，我没有发现您的回答很有帮助，IrrationalPerson。 @knedlsepp：我正在使用 in 来解决热泵中的热交换器。我有点在其他求解器中使用它。我猜这个函数至少被使用了 30.000 次，很可能更多（我正在做一年的每小时计算）。如前所述，x、b、c 和 d 是已知的。
您能否编辑您的帖子，为x,b,c,d 提供一些现实价值？我的第一个天真的想法是尝试使用fzero() 来查找(a-b-c+d)/ln((a-b)/(c-d)) - x 的根，但这可能会导致我尝试过的一些随机参数值的答案很复杂。
@ROLF：您是否已经尝试过评估lambertw 的向量而不是单个值？您可以使用lambertw(1:100) 代替for i = 1:100, lambertw(i), end，这比原来的速度快了大约30 倍。对于 30.000 个值，这大约需要 30 秒。
添加另一个想法：lambertw 在 matlab 中使用符号数学，这是一个巨大的开销。切换到数字实现。 octave 版本可能是最容易移植的：octave-specfun.sourcearchive.com/documentation/1.0.9-1/… 没有对代码进行基准测试，但只有 10 次迭代且没有其他循环的 for 循环应该很快。

标签： matlab logarithm wolframalpha natural-logarithm

【解决方案1】：

另一个选项基于更简单的Wright ω function：

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

假设d ~= c + x.*wrightOmega(log(-(c-d)./x) - (c-d)./x)（即，d ~= c+b-a、x 在这种情况下是0/0）。这相当于Lambert W function的主分支，W₀，我认为是你想要的解决方案分支。

与lambertW 一样，符号数学工具箱中也有一个wrightOmega 函数。不幸的是，对于大量输入，这可能也会很慢。但是，您可以使用我在 GitHub 上的 wrightOmegaq 进行复值浮点（双精度或单精度）输入。该函数更准确、完全矢量化，并且比使用内置 wrightOmega 进行浮点输入快三到四个数量级。

对于那些感兴趣的人，wrightOmegaq 基于这篇优秀的论文：

Piers W. Lawrence、Robert M. Corless 和 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.

该算法超越了 Cleve Moler 的 Lambert_W 中使用的 Halley's method 的三次收敛，并使用四阶收敛的求根方法（Fritsch, Shafer, & Crowley, 1973）收敛不超过两次迭代。

另外，为了进一步加快 Moler 的 Lambert_W 使用系列扩展，请参阅 my answer at Math.StackExchange。

【讨论】：

我对函数进行了一些测试以获得更好的概览：在区间 [-1/e,0] 中，您的起始值可以加快进程，因为 Moler 的实现需要大约6 次迭代（而不是两次）。在实践中，虽然它已经相当快了，因为它是完全矢量化的。你的回答有一点：不知道你说的FSC算法是不是更高阶，哈雷法是三阶法。
@knedlsepp：哈雷的方法是二阶Householder class method。收敛速度比阶数高一，因此哈雷方法具有三次收敛。 Lambert_W 函数将需要更多迭代才能达到某些临界值的给定容差。我的 Math.StackExchange 答案中的示例代码只是使用级数扩展的说明，而不是用于将速度与另一个函数进行比较的代码。该图显示更天真的方法使用更多迭代，尤其是在 -1/exp(1) 附近。
哦。 二阶Householder类方法：哈雷法？这可能会令人困惑，因为通常 nth-order 方法是关于收敛顺序的。（此外，FSC 论文并没有真正直接使用适用于这个问题的 Householder 方法，而是分别提供了三阶和四阶（收敛）的两种算法。）
@knedlsepp：我同意这件事令人困惑的顺序，但这不是我学到的。如果我想指定收敛，那么我总是明确地说 收敛速度 或 收敛顺序 是清楚的。而且，正如我的回答所说，它是劳伦斯等人。引用 FSC 论文作为我在 wrightOmegaq 中实现的迭代算法的基础的论文。
@ROLF：您可以忽略上面的 cmets - 它们与您的实际问题无关，仅与我回答中的措辞有关。从您的问题中我可以看出，Wright omega 函数应该可以正常工作。我尝试了一些数值并在顶部修改了我的公式。您应该插入一些您自己的值，并使用lambertw 将我上面的公式与您问题中的公式进行比较。两者应该返回非常相似的结果。 Matlab 的wrightOmega 可能会由于数值错误而返回一个很小的虚部（使用real 来消除它）。 wrightOmegaq 没有这个问题。

【解决方案2】：

两个（可组合的）选项：

您的脚本是否已经矢量化？为多个参数评估函数。执行for i = 1:100, a(i)=lambertw(rhs(i)); end 比a=lambertw(rhs) 慢。
如果您正在处理 LambertW 的实值分支（即您的参数在区间 [-1/e, inf) 中），您可以使用由 Cleve Moler 在 @ 上提交的 Lambert_W 的实现987654321@.

【讨论】：

Cleve Moler 的 Lambert_W 函数很好。我将推荐使用 exp 和 fzero 的组合，但 Moler 的做法基本相同，但使用的是三次方方法，因此它可能更快。
@TroyHaskin：是的，快速浏览了他关于 mathworks 的博客文章。不过，我看到使用比牛顿迭代更高阶的方法是很不寻常的。
@knedlsepp：感谢您的回复。请查看我更新的问题。
@ROLF：我不知道你更新后的问题应该告诉我什么？

【解决方案3】：

您知道每个时间步长热交换器两侧的质量流量吗？如果是，则温度“a”可以通过不需要任何迭代的“有效性-NTU”方法而不是 LMTD 方法来解决。参考：例如http://ceng.tu.edu.iq/ched/images/lectures/chem-lec/st3/c2/Lec23.pdf

【讨论】：