newton fractal: mathoverflow error

Question

I have written a python code for a newton method in 1D and want to use it to calculate the newton fractal for the function

The basic python code I am using is this:

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()

My newton()-function returns the number of iterations it needed to find an approximation xk such that |f(xk)|<1e-10.

I tested this code with and it worked, but when I use the actual function I would like to use, namely , I get an Overflow Error, "OverflowError: math range error". This is my code for the function f:

import cmath as cm

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

I don't really know how to debug this. I tried to convert to double precision, but my web research suggested that python is already using double precision. That's about the only idea I have for solving this.

Does anyone have an idea for what to do here? Thanks!

Edit:

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

Obviously the problem is with cmath.cos(). If you look at x when it fails, you see that it is

(8996.29520764-8256.38535552j)

Which is rather large.

---

Aside: You say you're not sure how to debug. There are lots of ways, but one easy way is to surround the statement that fails with a try...except block and then examine the value of x when an exception is encountered - e.g.

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

---

Using tigonometric substitutions for the complex maths cos function doesn't help either:

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

fails with a similar error, as does:

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

The latter fails because the hyperbolic functions fail. When you think about it - it is not very surprising that exp(8000) or so fails - that's a huge number. The largest that a 32 bit double can hold is about math.exp(709).

How does x get so large?

Your problem is that grad(x) is very small at some points in your function, causing xk to explode in value. This hapens as xk tends towards (0 + 0j) - grad tends to something very large!

You can see this by inserting print statements in the loop where you update xk

I'm not sure off the top of my head of the maths around how you're going to be able to control this, but I don't think that simply altering the precision will help. You might consider looking at the behaviour of the function around its roots.

Your other function (the simple quartic one) does not exhibit this problematic behaviour.