牛顿分形:mathoverflow错误

Question

我在1D中为newton方法编写了一个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)| <1e-10.

我测试了这段代码 它工作,但当我使用我想使用的实际功能,即 ,我得到一个溢出错误,"OverflowError:数学范围错误".这是函数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

Answer 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函数也没有帮助：

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趋向于非常大的东西！

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

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

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