这是我为使用牛顿法制作分形而编写的一个小脚本。
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/