logo标签列表

Matlab：解决对数方程

matlablogarithmwolframalphanatural-logarithm
7

我有以下方程式需要解决，其中要涉及到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页面上更容易地看到它。

我使用了Matlab内置的lambertW函数来解决这个方程。这相当慢，并且是我的脚本中的瓶颈。有没有另一种更快的方法来做到这一点？它不必精确到第10位小数。

编辑： 我不知道这个方程式有多难求。这里是一个说明我的问题的图片。温度b-d加上LMTD在每个时间步长中都会变化，但是是已知的。热从红线（CO2）传递到蓝线（水）。我需要找到温度“a”。我不知道这是如此难以计算！：P enter image description here

你需要多经常以及在什么情境下解决这个方程来求 a？是作为其他某个求解器的一部分吗？你想要评估的 x，b，c，... 的值是否已知？ - knedlsepp
我对 lambert W 函数一无所知。不过，IrrationalPerson，我觉得你的回答并不是很有帮助。@knedlsepp：我正在使用它来解决热泵中的热交换器问题。我在某个其他求解器中也使用了它。我想这个函数至少被使用了 30,000 次，很可能更多（我正在为一年的每小时进行计算）。如提到的，x、b、c 和 d 都是已知的。 - ROLF
你能否编辑你的帖子，提供一些x，b，c，d的实际值？我的第一个天真想法是尝试使用fzero()来找到(a-b-c+d)/ln((a-b)/(c-d)) - x的根，但对于我尝试过的一些随机参数值，这可能会导致复杂的答案。 - eigenchris
1
@ROLF：你是否已经尝试过对向量而非单个值进行lambertw的评估？你可以使用lambertw(1:100)代替for i = 1:100, lambertw(i), end，这样速度大约快了30倍。对于30000个值，这只需要大约30秒。 - knedlsepp
添加另一个想法：在Matlab中使用符号数学的lambertw是一个巨大的开销。切换到数字实现。Octave版本可能是最容易移植的：http://octave-specfun.sourcearchive.com/documentation/1.0.9-1/lambertw_8m-source.html 没有对代码进行基准测试，但仅使用10次迭代和没有其他循环的for循环应该很快。 - Daniel
显示剩余4条评论
3个回答
4

另一种选择基于更简单的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函数的主分支W0，我认为这是您想要的解决方案分支。

与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方法的立方收敛，并使用具有四阶收敛性 (Fritsch、Shafer和Crowley，1973) 的根查找方法，在不超过两次迭代中收敛。

此外，为了进一步加速Moler的 Lambert_W，可以使用级数展开，请参见 我在Math.StackExchange上的回答。

我对函数进行了一些测试，以便更好地了解情况：在区间[-1/e,0]中，您的起始值可以加快进程，因为Moler的实现需要约6次迭代（而不是两次）。但在实践中，它已经非常快速，因为它是完全向量化的。关于您的答案有一件事：我不知道您提到的FSC算法是否具有更高的阶数，但Halley方法是三阶方法。 - knedlsepp
1
@knedlsepp：哈雷方法是一个二阶豪斯霍尔德类方法。收敛速度比阶数高一级，因此哈雷方法具有立方收敛性。对于某些关键值，Lambert_W函数需要更多的迭代才能达到给定的容差。我在我的Math.StackExchange答案中提供的示例代码仅是使用级数展开的说明，而不是用于与另一个函数比较速度的代码。图表显示，更为朴素的方法使用更多的迭代，特别是在-1/exp(1)附近。 - horchler
@horchler：我对数学不是很感兴趣，所以我不理解讨论或您答案背后的数学。如果wrightOmega和wrightOmegaq都可以替换lambertw，那么我的理解是正确的吗？在我的代码中，我是否可以简单地用这两个其他函数之一替换lambertw（我知道wrightOmegaq不是内置的）？ - ROLF
1
@ROLF：你可以忽略上面的评论 - 它们与你实际问题无关，只是我回答中的措辞。从你的问题中最好的解释是Wright omega函数应该可以正常工作。我尝试了一些数值并稍微修改了我的公式。你应该插入一些自己的值，并使用lambertw将我的公式与你问题中的公式进行比较。两者应该返回非常相似的结果。Matlab的wrightOmega可能会由于数字误差而返回微小的虚部（使用real来消除它）。wrightOmegaq没有这个问题。 - horchler
谢谢！ :) 我不需要超过一到两位数字的精度。您上面的公式能否更改以提高准确性？ - ROLF
显示剩余12条评论
1

两个（可组合）选项：

  • 你的脚本已经向量化了吗？对多个参数求解函数值。执行for i=1：100，a(i)=lambertw(rhs(i)); end比a=lambertw(rhs)慢。
  • 如果你正在处理LambertW的实值分支（即你的参数在区间[-1/e, inf)），你可以使用由Cleve Moler提交的Lambert_W实现，在File Exchange上。
1
Cleve Moler 的 Lambert_W 函数很棒。我原本想推荐使用 exp 和 fzero 的组合，但是 Moler 实际上正在使用一个立方方法来完成基本相同的事情，所以它很可能更快。 - TroyHaskin
@TroyHaskin：是的，我快速浏览了他在mathworks上的博客文章。不过，对我来说看到比牛顿迭代更高阶的方法使用还是很不寻常的。 - knedlsepp
@knedlsepp：感谢您的回复。请查看我的更新问题。 - ROLF
@ROLF：我不知道你的更新问题应该告诉我什么？ - knedlsepp
0
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接