java - 尽管在早期迭代中运行，但程序似乎卡住了

Question

我正在编写一个程序，通过找到超稳定值，然后使用这些超稳定值的比率来计算常数，使用物流方程来计算费根鲍姆常数。

我对几乎所有值都使用 BigDecimals，以便在计算常量期间保持必要的精度水平。

我正在根据以下文件第 30-35 页上的 C++ 代码改编我的代码:http://webcache.googleusercontent.com/search?q=cache:xabTioRiF0IJ:home.simula.no/~logg/pub/reports/chaos_hw1.ps.gz+&cd=21&hl=en&ct=clnk&gl=us

我怀疑该程序的作用对我的问题是否重要。我运行了该程序，它似乎正在工作。我得到的前 4 个超稳定值和前 2 个 d 的输出是预期的，但在显示这 4 行后，程序似乎停止了。我没有遇到异常，但即使等待 30 分钟也没有输出更多计算结果。我无法弄清楚到底是什么导致了它，因为每行的计算时间应该大致相同，但显然不是。这是我的输出:

Feigenbaum constant calculation (using superstable points):
j       a           d
-----------------------------------------------------
1       2.0         N/A
2   3.23606797749979        N/A
4   3.4985616993277016  4.708943013540503
8   3.554640862768825   4.680770998010695

这是我的代码:

import java.math.*;


// If there is a stable cycle, the iterates of 1/2 converge to the cycle. 
// This was proved by Fatou and Julia. 
// (What's special about x = 1/2 is that it is the critical point, the point at which the logistic map's derivative is 0.)
// Source: http://classes.yale.edu/fractals/chaos/Cycles/LogisticCycles/CycleGeneology.html

public class Feigenbaum4
{
public static BigDecimal r[] = new BigDecimal[19];
public static int iter = 0;
public static int iter1 = 20; // Iterations for tolerance level 1
public static int iter2 = 10; // Iterations for tolerance level 2
public static BigDecimal tol1 = new BigDecimal("2E-31"); // Tolerance for convergence level 1
public static BigDecimal tol2 = new BigDecimal("2E-27"); // Tolerance for convergence level 2
public static BigDecimal step = new BigDecimal("0.01"); // step when looking for second superstable a
public static BigDecimal x0 = new BigDecimal(".5");
public static BigDecimal aZero = new BigDecimal("2.0");

public static void main(String [] args)
{
    System.out.println("Feigenbaum constant calculation (using superstable points):");
    System.out.println("j\t\ta\t\t\td");
    System.out.println("-----------------------------------------------------");

    int n = 20;
    if (FindFirstTwo())
    {
        FindRoots(n);
    }

}

public static BigDecimal F(BigDecimal a, BigDecimal x)
{
    BigDecimal temp = new BigDecimal("1");
    temp = temp.subtract(x);
    BigDecimal ans = (a.multiply(x.multiply(temp)));
    return ans;
}

public static BigDecimal Dfdx(BigDecimal a, BigDecimal x)
{
    BigDecimal ans = (a.subtract(x.multiply(a.multiply(new BigDecimal("2")))));
    return ans;
}

public static BigDecimal Dfda(BigDecimal x)
{
    BigDecimal temp = new BigDecimal("1");
    temp = temp.subtract(x);
    BigDecimal ans = (x.multiply(temp));
    return ans;
}

public static BigDecimal NewtonStep(BigDecimal a, BigDecimal x, int n)
{
    // This function returns the Newton step for finding the root, a,
    // of fn(x,a) - x = 0 for a fixed x = X

    BigDecimal fval = F(a, x);
    BigDecimal dval = Dfda(x);

    for (int i = 1; i < n; i++)
    {
        dval = Dfda(fval).add(Dfdx(a, fval).multiply(dval));
        fval = F(a, fval);
    }

    BigDecimal ans = fval.subtract(x);
    ans = ans.divide(dval, MathContext.DECIMAL64);
    ans = ans.negate();
    return ans;

}

public static BigDecimal Root(BigDecimal a0, int n)
{
    // Find the root a of fn(x,a) - x = 0 for fixed x = X
    // with Newton’s method. The initial guess is a0.
    //
    // On return iter is the number of iterations if
    // the root was found. If not, iter is -1.

    BigDecimal a = a0;
    BigDecimal a_old = a0;
    BigDecimal ans;

    // First iter1 iterations with a stricter criterion,
    // tol1 < tol2

    for (iter = 0; iter < iter1; iter++)
    {
        a = a.add(NewtonStep(a, x0, n));

        // check for convergence
        BigDecimal temp = a.subtract(a_old);
        temp = temp.divide(a_old, MathContext.DECIMAL64);
        ans = temp.abs();

        if (ans.compareTo(tol1) < 0)
        {
            return a;
        }

        a_old = a;
    }

    // If this doesn't work, do another iter2 iterations 
    // with the larger tolerance tol2
    for (; iter < (iter1 + iter2); iter++)
    {
        a = a.add(NewtonStep(a, x0, n));

        // check for convergence
        BigDecimal temp = a.subtract(a_old);
        temp = temp.divide(a_old, MathContext.DECIMAL64);
        ans = temp.abs();

        if (ans.compareTo(tol2) < 0)
        {
            return a;
        }

        a_old = a;
    }

    BigDecimal temp2 = a.subtract(a_old);
    temp2 = temp2.divide(a_old, MathContext.DECIMAL64);
    ans = temp2.abs();

    // If not out at this point, iterations did not converge
    System.out.println("Error: Iterations did not converge,");
    System.out.println("residual = " + ans.toString());

    iter = -1;

    return a;
}

public static boolean FindFirstTwo()
{
    BigDecimal guess = aZero;
    BigDecimal r0;
    BigDecimal r1;

    while (true)
    {
        r0 = Root(guess, 1);
        r1 = Root(guess, 2);

        if (iter == -1)
        {
            System.out.println("Error: Unable to find first two superstable orbits");
            return false;
        }

        BigDecimal temp = r0.add(tol1.multiply(new BigDecimal ("2")));
        if (temp.compareTo(r1) < 0)
        {
            System.out.println("1\t\t" + r0.doubleValue() + "\t\t\tN/A");
            System.out.println("2\t" + r1.doubleValue() + "\t\tN/A");

            r[0] = r0;
            r[1] = r1;

            return true;
        }

        guess = guess.add(step);

    }


}

public static void FindRoots(int n)
{
    int n1 = 4;
    BigDecimal delta = new BigDecimal(4.0);
    BigDecimal guess;

    for (int i = 2; i < n; i++)
    {
        // Computation

        BigDecimal temp = (r[i-1].subtract(r[i-2])).divide(delta, MathContext.DECIMAL64);
        guess = r[i-1].add(temp);
        r[i] = Root(guess, n1);
        BigDecimal temp2 = r[i-1].subtract(r[i-2]);
        BigDecimal temp3 = r[i].subtract(r[i-1]);
        delta = temp2.divide(temp3, MathContext.DECIMAL64);

        // Output

        System.out.println(n1 + "\t" + r[i].doubleValue() + "\t" + delta.doubleValue());

        // Step to next superstable orbit

        n1 = n1 * 2;
    }
}

}

编辑:Phil Steitz 的回答基本上解决了我的问题。我查看了一些线程转储，并做了一些研究来尝试理解它们，并使用调试信息编译了我的程序，我发现主线程停顿在以下行:

dval = Dfda(fval).add(Dfdx(a, fval).multiply(dval));

正如 Phil Steit 所说，通过使用

MathContext.DECIMAL128

不仅在这一行:

 dval = Dfda(fval).add(Dfdx(a, fval).multiply(dval));

而且在方法 F、Dfda 和 Dfdx 中的乘法运算中，我也能够让我的代码正常工作。

我使用 DECIMAL128 是因为较小的精度使计算不起作用，因为我将它们与如此低的数字进行比较以进行容差检查。

Answer 1

我认为这里发生的情况是，当 n 大于 10 时，您的 NewtonStep 方法会变得非常慢，因为您的任何 multiply 调用都不会限制规模通过提供 MathContext。当未提供 MathContext 时，乘法的结果将获得被乘数小数位的总和。使用上面的代码，NewtonStep 中 for 循环内的 dval 和 fval 对于较大的 n 来说会变得非常大，导致该方法以及以下方法中的乘法非常慢它呼唤。尝试在乘法激活中指定 MathContext.DECIMAL64 (或其他内容)，就像对除法所做的那样。