牛顿分形：数学溢出错误答案

【问题标题】：newton fractal: mathoverflow error牛顿分形：数学溢出错误
【发布时间】：2015-04-20 15:43:07
【问题描述】：

我已经为一维牛顿法编写了一个 python 代码，并想用它来计算函数的牛顿分形

我使用的基本python代码是这样的：

error = 1e-10
resolution = 100

x_range = np.linspace(-2,2,resolution)
y_range = np.linspace(-1,1,resolution)

fraktal = np.zeros(shape=(resolution,resolution))


for i in range(resolution):
    for j in range(resolution):
        x = x_range[i]
        y = y_range[j]
        z = complex(x,y)
        fraktal[i,j] = newton(z,error)

plt.imshow(fraktal)
plt.show()

我的 newton() 函数返回它需要的迭代次数来找到一个近似值 xk 使得 |f(xk)|

我用测试了这段代码，它可以工作，但是当我使用我想使用的实际函数时，即，我得到一个溢出错误，“溢出错误：数学范围错误”。这是我的函数 f 的代码：

import cmath as cm

def f(x):
    y = pow(x,4)*cm.cos(x) - 1
    return y

我真的不知道如何调试它。我试图转换为双精度，但我的网络研究表明 python 已经在使用双精度。这是我解决这个问题的唯一想法。

有人知道在这里做什么吗？谢谢！

编辑：

def grad(x):
    h = 1e-6
    y = (f(x+h)-f(x-h))/(2*h)
    return y

def newton(x0,error):
    k = 1

    xk = x0

    while 1:
        xk = xk - f(xk)/grad(xk)

        err = abs(f(xk))

        if err < error:
            return k
            break
        if k > 100:
            return 100
            break

        k = k+1

【问题讨论】：

我的朋友这里没有乳胶。不过使用CodeCogs 会很好看。
分享newton功能码怎么样？
好的，我编辑了 CodeCogs 并添加了牛顿方法的代码
我没有适合您的解决方案，但在我看来，由于迭代算法中评估的值很大，您会遇到此错误。请参阅这个相关问题：stackoverflow.com/questions/4050907/… 其中 exp 函数是罪魁祸首（对你来说，由于那些复杂的余弦，它可能是相同的）。每当我使用exp 时，我都会尝试检查指数是否太大（大于 709.78）

标签： python python-2.7 numpy

【解决方案1】：

显然问题出在cmath.cos()。如果你在失败时查看x，你会发现它是

(8996.29520764-8256.38535552j)

相当大。

旁白：你说你不确定如何调试。有很多方法，但一种简单的方法是用 try...except 块包围失败的语句，然后在遇到异常时检查 x 的值 - 例如

def f(x):
    try:
        y = pow(x,4)*cm.cos(x) - 1
    except OverflowError:
        print x
    return y

对复杂的数学 cos 函数使用 tigonometric 替换也无济于事：

j = cm.sqrt(-1)
cos_part = (cm.exp(j * x) + cm.exp(-1*j*x))/2

失败并出现类似的错误，如下所示：

cos_part = math.cos(x.real) * math.cosh(x.imag) - j * math.sin(x.real) * math.sinh(x.imag)

后者失败是因为双曲函数失败。仔细想想 - exp(8000) 左右失败并不奇怪 - 这是一个巨大数字。一个 32 位的double 可以容纳的最大值约为math.exp(709)。

x 怎么变得这么大？

您的问题是 grad(x) 在您的函数中的某些点非常小，导致 xk 价值爆炸。这是因为xk 趋向于(0 + 0j) - grad 趋向于非常大的东西！

您可以通过在更新 xk 的循环中插入 print 语句来看到这一点

我不确定您将如何控制这一点，但我认为简单地改变精度不会有帮助。您可以考虑查看函数围绕其根的行为。

您的其他函数（简单的四次函数）没有表现出这种有问题的行为。

【讨论】：