java - 在 Java 中从 Matlab 计算复杂的算术表达式

Question

我正在使用 Matlab 推导一些符号方程(带变量)。这些是非常困难的函数，包括 log、tanh、atanh 函数:

import static java.lang.Math.log;
import static java.lang.Math.tanh;

public class MathTest {

    private static double coeff(double accel_Vmax, double accel_amin, double accel_amax, double accel_x0, double accel_v0) {
        double accel_Vmax2 = Math.pow(accel_Vmax, 2);
        return -(1.6658627956859751925e-24*(618147116152812568000.0*accel_Vmax2*log(1.0 - 1.0*tanh((0.071428571428571428571*(660.0*accel_amax - 165.0*accel_amin + 14.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) + 3.3331462164783108864e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.14285714285714285714*(180.0*accel_amax - 45.0*accel_amin + 7.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) - 6.0723863811437384517e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.071428571428571428571*(180.0*accel_amax - 45.0*accel_amin + 14.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) - 3.9876549298255135726e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.14285714285714285714*(360.0*accel_amax - 90.0*accel_amin + 7.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) - 1.8180797558373403944e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.5*(60.0*accel_amax - 15.0*accel_amin + 2.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) - 6.002895332014512297e23*accel_x0*accel_amin + 1.2005790664029024594e24*accel_x0*accel_amax + 262429425.0*accel_Vmax*accel_amin + 2.4362268723458870243e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.071428571428571428571*(780.0*accel_amax - 195.0*accel_amin + 14.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) + 8.8358676058937131753e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.14285714285714285714*(60.0*accel_amax - 15.0*accel_amin + 7.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) - 4.2543066242278358162e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.071428571428571428571*(60.0*accel_amax - 15.0*accel_amin + 14.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) - 3.1634587719706994218e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.14285714285714285714*(120.0*accel_amax - 30.0*accel_amin + 7.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) - 3.6361595097781225383e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.14285714285714285714*(240.0*accel_amax - 60.0*accel_amin + 7.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) - 472700736468785811844.0*accel_Vmax2*log(1.0 - 1.0*tanh((60.0*accel_amax - 15.0*accel_amin + accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax)))/accel_Vmax)) + 735312255896455259166.0*accel_Vmax2*log(-(4.0*(accel_v0*accel_amin - 2.0*accel_v0*accel_amax - 1.0*accel_Vmax*accel_amin + 2.0*accel_Vmax*accel_amax))/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax)) + 101004431699873120300.0*accel_Vmax2*log(1.0 - 1.0*tanh((0.071428571428571428571*(540.0*accel_amax - 135.0*accel_amin + 14.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) + 3.5270747250077918894e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.14285714285714285714*(300.0*accel_amax - 75.0*accel_amin + 7.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax)) + 3.817967485777192084e21*accel_Vmax2*log(1.0 - 1.0*tanh((0.071428571428571428571*(300.0*accel_amax - 75.0*accel_amin + 14.0*accel_Vmax*atanh((4.0*accel_v0*accel_amin - 8.0*accel_v0*accel_amax - 3.0*accel_Vmax*accel_amin + 4.0*accel_Vmax*accel_amax)/(accel_Vmax*accel_amin - 4.0*accel_Vmax*accel_amax))))/accel_Vmax))))/(accel_amin - 2.0*accel_amax);
    }

    private static double atanh(double x) {
        return 0.5 * log((1.0 + x) / (1.0 - x));
    }

    public static void main(String[] args) {
        System.out.println("" + coeff(40, 2, 15, 0,0));
    }
}

如您所见，coeff() 函数中的表达式非常糟糕。存在 double 变点算术失败的重复模式:

log(1.0 - 1.0*tanh(...))

通常 tanh(...) 的结果类似于 0.9998232、0.9892348、0.99999999843、0.999999999999... 但总是在某些时候它输出 1 而 log() 变为 NaN。

我想知道是否有一种解决方案可以在不从Java调用Matlab API的情况下计算这个表达式。我的目标是性能，无法在生产中使用 Matlab。它是可以实现的还是超出范围的？

Answer 1

您可以将表达式 1-tanh(x) 简化为 2/(e^2x - 1)。这应该可以解决对数内部发生的灾难性抵消。 (您仍然需要担心 x 接近 0 处的零分母，但简单的特殊情况就足够了。)

编辑:更进一步，您可以利用 Math.log1p()，它是为此类事情量身定制的:

  log(1 - tanh(x))
= log(1 - (e^2x - 1)/(e^2x + 1))
= log(2/(e^2x + 1))
= log(2) - log(e^2x + 1)
= log(2) - log1p(e^2x)

看起来你的方程比 log(1-tanh(x)) 稍微复杂一点，但基本方法是，尽可能简化，然后利用数学函数，例如Math.log1p() 以避免不稳定区域。 (一般来说，您希望数字“接近”的唯一数字是零。)