朗之万函数发生在与弹性体和顺磁性材料相关的计算中。只要解决了原点附近的减法抵消问题,例如通过在该区域中使用极小极大多项式近似,就可以轻松准确且高性能地计算它。
为了完全消除对消问题,选择一个良好的多项式近似转换点很重要;一些实现统一切换,这是次优的。更接近的切换点是可取的。忠实地实现指数函数,通常由现代数学库提供,最大误差为ulp 可以直接实现。示例性 C 实现如下所示。
然而,许多应用程序需要反朗之万函数,没有封闭形式的表示;需要数值近似。文献中使用的大多数近似值都非常不准确,相对误差在 0.05% 到 2% 的范围内。有用的概述由以下人员提供:
Benjamin C. Marchi 和 Ellen M. Arruda,“逆朗之万近似的误差最小化方法”。流变学报(2015) 54:887-902
Martin Kröger,“朗之万和布里渊逆函数的简单、可接受和准确的近似值,与强聚合物变形和流动相关。” 非牛顿流体力学杂志223 (2015) 77-87
最近的工作表明,反朗之万函数的计算缺乏准确性对其使用的科学计算的整体准确性产生负面影响:
Amine Ammar,“逆朗之万近似对非线性稀聚合物福克-普朗克方程解的影响。” 非牛顿流体力学杂志231(2016):1-5。
因此,非常需要准确计算反朗之万函数。除了准确性之外,一些用例还受益于特定的功能形式,以便于替换、区分和集成。抛开这些额外的要求,专注于最小化相对误差,反朗之万函数如何计算准确(接近机器精度)并具有良好的性能?
#include <math.h> // fabs, exp, copysign
/* compute Langevin function L(a) = coth(a) - 1/a. Max error typically < 1.3 ulp */
float langevinf (float a)
{
float r, t;
t = fabsf (a);
if (t > 1.8125f) {
float e;
e = 1.0f / (expf (2.0f * t) - 1.0f);
r = (1.0f - 1.0f / t) + 2.0f * e;
r = copysignf (r, a);
} else {
float s;
s = a * a;
r = 7.70960469e-8f;
r = r * s - 1.65101926e-6f;
r = r * s + 2.03457112e-5f;
r = r * s - 2.10521728e-4f;
r = r * s + 2.11580913e-3f;
r = r * s - 2.22220998e-2f;
r = r * s + 8.33333284e-2f;
r = r * a + 0.25f * a;
}
return r;
}
/* compute Langevin function L(a) = coth(a) - 1/a. Max error typically < 1.3 ulp */
double langevin (double a)
{
double r, t;
t = fabs (a);
if (t > 1.8125) {
double e;
e = 1.0 / (exp (2.0 * t) - 1.0);
r = (1.0 - 1.0 / t) + 2.0 * e;
r = copysign (r, a);
} else {
double s, t, u;
s = a * a;
t = s * s;
r = 2.6224020345921745e-16 * s - 9.0382507125284313e-15;
u = 1.6058476918889783e-13 * s - 2.0544032928760129e-12;
r = r * t + u;
u = 2.2318295632954910e-11 * s - 2.2645716105162602e-10;
r = r * t + u;
u = 2.2484176131435534e-09 * s - 2.2212056151537571e-08;
r = r * t + u;
u = 2.1925751794205211e-07 * s - 2.1644032423757088e-06;
r = r * t + u;
u = 2.1377798795067705e-05 * s - 2.1164021156459166e-04;
r = r * t + u;
u = 2.1164021163937877e-03 * s - 2.2222222222221859e-02;
r = r * t + u;
r = r * s + 8.3333333333333329e-02;
r = r * a + 0.25 * a;
}
return r;
}