Mandelbrot Set in R

r blogger の ”A 3D Version of R’s curve() Function” という投稿のコメント欄にある lattice ライブラリの例にインスパイアされて、マンデルブロー集合(Mandelbrot Set)を R & ggplot2 /lattice で 2D/3D グラフを書いてみた。

Definition
Mandelbort 集合z_{n+1} = z_n^2 + c, c in C(複素数), z_0 = 0
で定義される漸化式で、n -> ∞ の時に収束するような c 全体の集合。

Program

library(ggplot2)
library(lattice)
g <- function(x0, y0) {
  x <- 0
  y <- 0
  for (i in 1:20) {
    xtemp <- x ^ 2 - y ^ 2 + x0
    y <- 2 * x * y + y0
    x <- xtemp
  }
  exp(-(x^2+y^2))
}
m <- 1200
X <- seq(-1.8, 0.6, length.out=m)
Y <- seq(-1.2, 1.2, length.out=m)
grid   <- expand.grid(x=X, y=Y)
grid$z <- g(grid$x, grid$y)

# lattice 3D plot
wireframe(z ~ x*y, grid, shade=TRUE)
# ggplot2 2D plot
ggplot(grid[! is.na(grid$z), ], aes(x=x, y=y, fill=z)) + geom_tile()

MEMO1
wikipedia にある Mandelbrot のコードとインスパイアされた r-blogger のコードを参考にさせてもらった。
なお wikipedia では漸化式を C^1 上で計算しているが、自分のコードは 実部、虚部の R^2 で計算。

MEMO2
wikipedia のコードでは次のような箇所がある。

X[,,k] <- exp(-abs(Z))

最初意味が分からなかったが、Z が発散するなら漸化式からその絶対値は ∞ になる(十分に大きくなる)から、絶対値の負の値を指数関数に与えることで、収束/発散を

  • 実質0 = 発散
  • 0 より大きい = 収束

と判定している(正しく理解しているなら)。

MEMO3
Mandelbrot Set の名前を知ってはや15年くらいは経過している。
MathematicaMaple で簡単に書けるのは昔から知っていたが、自分でプログラムで書いたのは初。

Advertisements
Tagged with: , , , , ,
Posted in algorithm, R
2 comments on “Mandelbrot Set in R
  1. […] 以前書いたコードをそのまま流用したのが以下 […]

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Archives
  • RT @__apf__: How to write a research paper: a guide for software engineers & practitioners. docs.google.com/presentation/d… /cc @inwyrd 6 months ago
  • RT @HayatoChiba: 昔、自然と対話しながら数学に打ち込んだら何かを悟れるのではと思いたち、専門書1つだけ持ってパワースポットで名高い奈良の山奥に1週間籠ったことがある。しかし泊まった民宿にドカベンが全巻揃っていたため、水島新司と対話しただけで1週間過ぎた。 それ… 6 months ago
  • RT @googlecloud: Ever wonder what underwater fiber optic internet cables look like? Look no further than this deep dive w/ @NatAndLo: https… 6 months ago
  • @ijin UTC+01:00 な時間帯で生活しています、、、 1 year ago
  • RT @mattcutts: Google's world-class Site Reliability Engineering team wrote a new book: amazon.com/Site-Reliabili… It's about managing produc… 1 year ago
%d bloggers like this: