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

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 linspace 的每个元素以应用 newton() 函数,该函数作用于单个猜测而不是整个数组。

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

P.S.:如果您想在不等待太长时间的情况下尝试代码,最好使用 100×100。

额外背景:
参见牛顿法寻找多项式的零点。
分形的基本思想是在复平面中测试猜测并计算迭代次数以收敛到零。这就是 newton() 中的递归,它最终返回步数。复平面中的猜测代表图片中的一个像素,由收敛的步数着色。从一个简单的算法中,你会得到这些漂亮的分形。

最佳答案

我使用了 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 时,这实际上是加速的一个非常重要的部分。

第二个加速来自仅对未收敛部分进行计算。劳里茨为此使用了掩码阵列,然后才意识到它们正在减慢速度。我已经有一段时间没有看过它们了,但我确实记得掩码数组在过去是缓慢的根源。我相信这是因为它们主要是在 numpy 数组上用纯 Python 实现的,而不是像 numpy 数组一样几乎完全用 C 编写。

关于python - 如何使用numpy数组加速分形生成?,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/17393592/

10-16 13:04
查看更多