7

这是我为使用牛顿法制作分形而编写的一个小脚本。

import numpy as np
import matplotlib.pyplot as plt

f = np.poly1d([1,0,0,-1]) # x^3 - 1
fp = np.polyder(f)

def newton(i, guess):
    if abs(f(guess)) > .00001:
        return newton(i+1, guess - f(guess)/fp(guess))
    else:
        return i

pic = []
for y in np.linspace(-10,10, 1000):
    pic.append( [newton(0,x+y*1j) for x in np.linspace(-10,10,1000)] )

plt.imshow(pic)
plt.show()

我正在使用 numpy 数组，但仍然循环遍历 1000×1000 linspaces 的每个元素以应用该newton()函数，该函数作用于单个猜测而不是整个数组。

我的问题是：如何改变我的方法以更好地利用 numpy 数组的优势？

PS：如果你想在不等待太久的情况下尝试代码，最好使用 100×100。

额外背景：
参见牛顿法找到多项式的零点。
分形的基本思想是在复平面中测试猜测并计算迭代次数以收敛到零。这就是 in 中的递归newton()，它最终返回步数。复平面中的猜测代表图片中的一个像素，按收敛的步数着色。通过一个简单的算法，您可以获得这些美丽的分形。

4

4 回答 4

5

我从 Lauritz V. Thaulow 的代码开始工作，并且能够通过以下代码获得相当显着的加速：

import numpy as np
import matplotlib.pyplot as plt
from itertools import count

def newton_fractal(xmin, xmax, ymin, ymax, xres, yres):
    yarr, xarr = np.meshgrid(np.linspace(xmin, xmax, xres), \
                             np.linspace(ymin, ymax, yres) * 1j)
    arr = yarr + xarr
    ydim, xdim = arr.shape
    arr = arr.flatten()
    f = np.poly1d([1,0,0,-1]) # x^3 - 1
    fp = np.polyder(f)
    counts = np.zeros(shape=arr.shape)
    unconverged = np.ones(shape=arr.shape, dtype=bool)
    indices = np.arange(len(arr))
    for i in count():
        f_g = f(arr[unconverged])
        new_unconverged = np.abs(f_g) > 0.00001
        counts[indices[unconverged][~new_unconverged]] = i
        if not np.any(new_unconverged):
            return counts.reshape((ydim, xdim))
        unconverged[unconverged] = new_unconverged
        arr[unconverged] -= f_g[new_unconverged] / fp(arr[unconverged])

N = 1000
pic = newton_fractal(-10, 10, -10, 10, N, N)

plt.imshow(pic)
plt.show()

对于 N=1000，我使用 Lauritz 的代码获得 11.1 秒的时间，使用此代码获得 1.7 秒的时间。

这里有两个主要的加速。首先，我使用 meshgrid 来加速输入值的 numpy 数组的创建。当 N=1000 时，这实际上是加速的一个非常重要的部分。

第二个加速来自仅对未收敛部分进行计算。Lauritz 曾为此使用掩码数组，然后才意识到它们正在减慢速度。我已经有一段时间没有看过它们了，但我确实记得蒙面数组是过去缓慢的根源。我相信这是因为它们在很大程度上是在纯 Python 中通过 numpy 数组实现的，而不是像 numpy 数组那样几乎完全用 C 语言编写。

于 2013-07-01T03:52:50.523 回答
3

这是我的尝试。它快了大约 16 倍。

import numpy as np
import matplotlib.pyplot as plt
from itertools import count

def newton_fractal(xmin, xmax, ymin, ymax, xres, yres):
    arr = np.array([[x + y * 1j for x in np.linspace(xmin, xmax, xres)] \
        for y in np.linspace(ymin, ymax, yres)], dtype="complex")
    f = np.poly1d([1,0,0,-1]) # x^3 - 1
    fp = np.polyder(f)
    counts = np.zeros(shape=arr.shape)
    for i in count():
        f_g = f(arr)
        converged = np.abs(f_g) <= 0.00001
        counts[np.where(np.logical_and(converged, counts == 0))] = i
        if np.all(converged):
            return counts
        arr -= f_g / fp(arr)

pic = newton_fractal(-10, 10, -10, 10, 100, 100)

plt.imshow(pic)
plt.show()

我不是一个 numpy 专家，我相信那些可以优化它的人，但在速度方面已经是一个巨大的改进。

编辑：原来屏蔽数组根本没有帮助，删除它们导致速度提高了 15% ，所以我从上述解决方案中删除了屏蔽数组。谁能解释为什么蒙面数组没有帮助？

于 2013-06-30T19:53:41.710 回答
3

我对牛顿函数进行了矢量化，得到了大约。200x200 点快 85 倍，500x500 点快 144 倍，1000x1000 点快 148 倍：

import numpy as np
import matplotlib.pyplot as plt

f = np.poly1d([1,0,0,-1]) # x^3 - 1
fp = np.polyder(f)
def newton(i, guess):            
    a = np.empty(guess.shape,dtype=int)
    a[:] = i
    j = np.abs(f(guess))>.00001
    if np.any(j):         
        a[j] = newton(i+1, guess[j] - np.divide(f(guess[j]),fp(guess[j])))        
    return a

npts = 1000
x = np.linspace(-10,10,npts)
y = np.linspace(-10,10,npts)
xx, yy = np.meshgrid(x, y)
pic = np.reshape(newton(0,np.ravel(xx+yy*1j)),[npts,npts])
plt.imshow(pic)
plt.show()
于 2013-07-01T02:53:23.167 回答
0

好的，我已经解决了 Justin Peel 代码中的无限循环，在代码中添加了最大迭代条件，现在代码绘制了 z^4-1 之类的多项式，并且它不会进入无限循环。如果有人知道如何改进此错误，请告诉我们。我的解决方案可能会减慢代码的执行速度，但它可以工作。这是代码：

    #!/usr/bin/python
    # -*- coding: utf-8 -*-

    import numpy as np
    import itertools
    import matplotlib.pyplot as plt

    __author__ = 'Tobal'
    __version__ = 1.0


    def newton_fractal(f, xmin, xmax, ymin, ymax, xres, yres, tolerance, maxiters):
        yarr, xarr = np.meshgrid(np.linspace(xmin, xmax, xres), np.linspace(ymin, ymax, yres) * 1j)
        arr = yarr + xarr
        ydim, xdim = arr.shape
        arr = arr.flatten()
        fp = np.polyder(f, m=1)
        counts = np.zeros(shape=arr.shape)
        unconverged = np.ones(shape=arr.shape, dtype=bool)
        indices = np.arange(len(arr))
        iters = 0
        for i in itertools.count():
            f_g = f(arr[unconverged])
            new_unconverged = np.abs(f_g) > tolerance
            counts[indices[unconverged][~new_unconverged]] = i
            if not np.any(new_unconverged) or iters >= maxiters:
                return counts.reshape((ydim, xdim))
            iters += 1
            unconverged[unconverged] = new_unconverged
            arr[unconverged] -= f_g[new_unconverged] / fp(arr[unconverged])


    pic = newton_fractal(np.poly1d([1., 0., 0., 0., -1.]), -10, 10, -10, 10, 1000, 1000, 0.00001, 1000)
    plt.imshow(pic, cmap=plt.get_cmap('gnuplot'))
    plt.title(r'$Newton^{\prime} s\;\;Fractal\;\;Of\;\;P\,(z) =z^4-1$', fontsize=18, y=1.03)
    plt.tight_layout()
    plt.show()

我正在将 Pycharm 5 与 Anaconda Python 3 一起使用，并且此 IDE 在代码中报告了一个警告，而 不是 np.any(new_unconverged)

预期的'type Union [ndarray，iterable]'，改为'bool'

这个警告也出现在 Justin Peel 的原始代码中。而且我不知道如何解决它。我对这个问题很感兴趣。 牛顿分形

于 2016-01-06T10:40:46.087 回答