2017-07-28 4 views
1

私はこれに関する他の質問を見てきましたが、どこが間違っているのか分かりません。目標は、-2≦x≦2かつ-2≦yの領域にまたがるN×Nグリッド上のc = x + iyのすべての値に対して反復を実行することによってマンデルブロ集合の画像を作るプログラムを書くことです≤2。マンデルブロ集合内の格子点が黒色で外側が白色である密度プロットを作成する。マンデルブロ集合内にあるとみなされる場合、z '= z^2 + cとして反復されたとき、zの大きさは決して2より大きくならない。Pythonで設定された単純なマンデルブロ

私のコードでは、ほとんど完全に黒い1つの小さなノッチがあります。将来の読者のために

from pylab import imshow,show,gray 
from numpy import zeros,linspace 

z = 0 + 0j 
n=100 

M = zeros([n,n],int) 
xvalues = linspace(-2,2,n) 
yvalues = linspace(-2,2,n) 

for x in xvalues: 
    for y in yvalues: 
     c = complex(x,y) 
     for i in range(100): 
      z = z*z + c 
      if abs(z) > 2.0: 
       M[y,x] = 1 
       break 

imshow(M,origin="lower") 
gray() 
show() 

、ここに私の新しいコードは見終わった方法です:

from pylab import imshow,show,gray 
from numpy import zeros,linspace 

n=1000 

M = zeros([n,n],int) 
xvalues = linspace(-2,2,n) 
yvalues = linspace(-2,2,n) 

for u,x in enumerate(xvalues): 
    for v,y in enumerate(yvalues): 
     z = 0 + 0j 
     c = complex(x,y) 
     for i in range(100): 
      z = z*z + c 
      if abs(z) > 2.0: 
       M[v,u] = 1 
       break 

imshow(M,origin="lower") 
gray() 
show() 

答えて

2

あなたのコードで問題のカップルがあります。

xvaluesyvaluesを使用してインデックスMをインデックスしていますが、これらのインデックスは0〜(n-1)の範囲のピクセルインデックス整数である必要があります。 Numpyはxvaluesyvaluesの浮動小数点数を整数に変換しますが、結果の数値は-2..2になるため、プロットされたピクセル数は少なくなります。画像は小さく、画像は小さくなります。双方向の負のインデックスはPythonで動作します。

必要なピクセルインデックスを取得する簡単な方法は、組み込みのPython関数enumerateを使用することですが、Numpy関数でこれを行うためにコードを再編成する方法があるかもしれません。

もう1つの問題は、ピクセルごとにzをゼロにリセットする必要があることです。現在、あなたのコードは、最後のzが前のピクセルのものであったものを再利用し、そのピクセルがマンデルブロセットにあった場合、zは大きすぎます。

修正されたバージョンのコードです。私はピラブを持っていないので、PILを使って簡単なビットマップビューアを書いた。私のshow機能の中でimg.save(filename)を呼び出すことによって、イメージをファイルに保存することができます。 PILは、ファイル名拡張子から正しいファイル形式を探し出します。

import numpy as np 
from PIL import Image 

def show(data): 
    img = Image.frombytes('1', data.shape[::-1], np.packbits(data, 1)) 
    img.show() 

n = 100 
maxiter = 100 

M = np.zeros([n, n], np.uint8) 
xvalues = np.linspace(-2, 2, n) 
yvalues = np.linspace(-2, 2, n) 

for u, x in enumerate(xvalues): 
    for v, y in enumerate(yvalues): 
     z = 0 
     c = complex(x, y) 
     for i in range(maxiter): 
      z = z*z + c 
      if abs(z) > 2.0: 
       M[v, u] = 1 
       break 

show(M) 

ここで、出力画像です:もちろん

B&W Mandelbrot


、あなたがnumpyの配列のインデックスを反復自分自身を見つけるたび、それは、あなたはそれが間違ってやっている兆候です。 Numpyを使用する主なポイントは、C言語で内部的に反復処理を行うことで、一度に配列全体の操作を実行できることです。上記のコードは、Numpy配列の代わりに単純なPythonリストを使用することもできます。

Numpyがほとんどのループを行うバージョンです。それはより多くのRAMを使用しますが、それは以前のバージョンよりおよそ2.5倍速く、やや短くなっています。

このコードでは、複数の2D配列を使用しています。 cにはすべての複合シード数が含まれており、コアマンデルブロ計算はzに実行されます。この計算はゼロに初期化されます。 maskは、マンデルブロ計算の実行場所を制御するブール値配列です。すべての項目は最初にTrueに設定され、各繰り返しで maskの項目は、の項目がマンデルブロセットからエスケープされた項目に対応し、Falseに設定されています。ポイントは、我々はそれが高価な平方根計算とabs関数呼び出しを回避するので、これは少し時間を節約する、というよりもabs(z) > 2.0z.real**2 + z.imag**2 > 4.0を使用して脱出したかどうかをテストするには

Mandelbrotセットをプロットするのに最終値maskを使用できますが、Mandelbrotのポイントを黒に設定するには、値を逆にする必要があります(1 - maskで行うことができます)。

import numpy as np 
from PIL import Image 

def show(data): 
    img = Image.frombytes('1', data.shape[::-1], np.packbits(data, 1)) 
    img.show() 
    img.save('mset.png') 

n = 100 
maxiter = 100 

a = np.linspace(-2, 2, n) 
c = a + 1.j * a[:, None] 
z = np.zeros((n, n), np.complex128) 
mask = np.ones((n, n), np.bool) 

for i in range(maxiter): 
    mask[mask] = z[mask].real**2 + z[mask].imag**2 < 4.0 
    z[mask] = z[mask]**2 + c[mask] 

show(1 - mask) 
+0

ありがとうございます!私はそれを働かせることができました。私はまだ列挙関数を学んでいないので、それは本当に便利です! –

+0

@NoraBaileyあなたが好きとうれしいです。 'enumerate'はプレーンなPythonで作業するときの非常に便利な関数ですが、Numpyには便利なことはまれです。なぜなら、通常、Numpyはあなたのためにループを見ています。 –

+0

@NoraBaileyあなたはこれが好きかもしれませんhttp://www.linuxtopia.org/online_books/programming_books/python_programming/python_ch20s03.htmlこれは、Pythonで他のリスト処理ユーティリティをカバーしています:) –

関連する問題