Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down on April 23, and its replacement, statalist.org is already up and running.


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: Mandelbrot sets


From   Scott Merryman <scott.merryman@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   st: Mandelbrot sets
Date   Wed, 14 Nov 2012 11:47:51 -0600

In the recent issue of The Stata News
(http://www.stata.com/stata-news/statanews.27.4.pdf) there is an
article by Hua Peng on generating and graphing Mandelbrot sets.

There does appear to be an error in the second block of code.  The
line res= J(201*201, 3, .) ,  I believe should be mandelbrot_set =
J(201*201, 3, .)

Warning - it can take a while for the graphs to be produced.

Scott


clear*
mata:
real scalar escape(complex scalar Z, C,
                   real scalar    R, maxiter)
{
real scalar i
for (i=0; i<maxiter; i++) {
   if (norm(Z) <= R) {
    Z = Z*Z+C
   }
else  return(i)
}
return(i)
}
end

mata:
mandelbrot_set = J(201*201, 3, .)
cnt = 1
C   = -0.8+0.156i
for(i=0; i<=200; i++) {
for(j=0; j<=200; j++) {
Z  = C(-2+i*4/200, -2+j*4/200)
n = escape(Z, C, 100, 400)
mandelbrot_set[cnt, .] = (i, j, n)
cnt++
}
}
end

getmata (x y escape)= mandelbrot_set
qui summ escape
local min = r(min)-1
local max = r(max)+1
tw contour escape y x, ccuts(`min'(1)`max') clegend(off)
graphr(m(zero)) xscale(off) xlab(,nogrid) yscale(off) ylab(,nogrid)
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index