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 列の linspace の各要素をループして、newton()配列全体ではなく単一の推測に作用する関数を適用します。

私の質問は次のとおりです。numpy 配列の利点をより有効に活用するために、アプローチをどのように変更できますか?

PS: あまり長く待たずにコードを試してみたい場合は、100 x 100 を使用することをお勧めします。

余分な背景:
多項式のゼロを見つけるためのニュートンの方法を参照してください。
フラクタルの基本的な考え方は、複素平面で推測をテストし、ゼロに収束するまでの反復回数を数えることです。それが の再帰でnewton()あり、最終的にステップ数を返します。複素平面の推定値は、画像内のピクセルを表し、収束までのステップ数で色分けされています。シンプルなアルゴリズムから、これらの美しいフラクタルが得られます。

4

4 に答える 4

5

私はLauritz V. Thoulowのコードから作業し、次のコードでかなりのスピードアップを得ることができました:

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 秒の時間が得られます。

ここでは主に 2 つの高速化が行われます。まず、meshgrid を使用して、numpy 配列の入力値の作成を高速化しました。これは、N=1000 の場合のスピードアップのかなり重要な部分です。

2 番目の高速化は、収束していない部分でのみ計算を行うことから得られます。ローリッツは、マスクされた配列が速度を低下させていることに気付く前に、これを使用していました。私はかなり長い間それらを見ていませんでしたが、過去にマスクされた配列が速度低下の原因だったことを覚えています。これは、numpy 配列のようにほぼ完全に C で記述されているのではなく、numpy 配列を介して純粋な Python で大部分が実装されているためだと思います。

于 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 に答える
0

OK、ジャスティン・ピールのコードで最大反復条件をコードに追加して無限ループを解決しました。コードは 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()

Anaconda Python 3 で Pycharm 5 を使用していますが、この IDE は np.any(new_unconverged) ではないコードで警告を報告します

'type Union[ndarray, iterable]' を予期していましたが、代わりに 'bool' を取得しました

この警告は、元の Justin Peel コードにも表示されます。そして、私はそれを解決する方法を知りません。私はこの問題に非常に興味があります。 ニュートンのフラクタル

于 2016-01-06T10:40:46.087 に答える