Sage Interactions - Fractal
goto interact main page
Contents
-
Sage Interactions - Fractal
- Mandelbrot's Fractal Binomial Distribution
- Fractals Generated By Digit Sets and Dilation Matrices
- Demonstrating that the Twin Dragon Matrix is likely to yield a Tiling of a Compact Interval of R^2 as k->infinity (It does!)
- Now in 3D
- Exploring Mandelbrot
- Mandelbrot & Julia Interact with variable exponent
- Sierpiński Triangle
Mandelbrot's Fractal Binomial Distribution
xxxxxxxxxxdef muk_plot(m0,k): """ Return a plot of the binomial fractal measure mu_k associated to m0, 1-m0, and k. """ k = int(k) m0 = float(m0) m1 = float(1 - m0) assert m0 > 0 and m1 > 0, "both must be positive" v = [(0,0)] t = 0 two = int(2) delta = float(1/2^k) multiplier = float(2^k) for i in [0..2^k-1]: t = i * delta phi1 = i.str(two).count("1") phi0 = k - phi1 y = m0^(phi0)*m1^(phi1)*multiplier v.append((t,y)) v.append((t+delta,y)) return vpretty_print(html("<h1>Mandelbrot's Fractal Binomial Measure</h1>"))def _(mu0=slider(0.0001,0.999,default=0.3), k=slider([1..14],default=3), thickness=slider([0.1,0.2,..,1.0],default=1.0)): v = muk_plot(mu0,k) line(v,thickness=thickness).show(xmin=0, xmax=1, ymin=0, figsize=[8,3])
Fractals Generated By Digit Sets and Dilation Matrices
(Sage Days 9 - Avra Laarakker)
Attempt at Generating all integer vectors with Digits D and Matrix A (How about vector([0,-1])?)
xxxxxxxxxxA = matrix([[1,1],[-1,1]])D = [vector([0,0]), vector([1,0])]def f(A = matrix([[1,1],[-1,1]]), D = '[[0,0],[1,0]]', k=[3..17]): print("Det = {}".format(A.det())) D = matrix(eval(D)).rows() def Dn(k): ans = [] for d in Tuples(D, k): s = sum(A^n*d[n] for n in range(k)) ans.append(s) return ans G = points([v.list() for v in Dn(k)],size=50) show(G, frame=True, axes=False)
Demonstrating that the Twin Dragon Matrix is likely to yield a Tiling of a Compact Interval of R^2 as k->infinity (It does!)
xxxxxxxxxxA = matrix([[1,1],[-1,1]])D = [vector([0,0]), vector([1,0])]def f(A = matrix([[1,1],[-1,1]]), D = '[[0,0],[1,0]]', k=[3..17]): print("Det = {}".format(A.det())) D = matrix(eval(D)).rows() def Dn(k): ans = [] for d in Tuples(D, k): s = sum(A^(-n)*d[n] for n in range(k)) ans.append(s) return ans G = points([v.list() for v in Dn(k)]) show(G, frame=True, axes=False)
Now in 3D
xxxxxxxxxxA = matrix([[0,0,2],[1,0,1],[0,1,-1]])D = '[[0,0,0],[1,0,0]]'def Dn(D,A,k): ans = [] for d in Tuples(D, k): s = sum(A^n*d[n] for n in range(k)) ans.append(s) return ansdef f(A = matrix([[0,0,2],[1,0,1],[0,1,-1]]), D = '[[0,0,0],[1,0,0]]', k=[3..15], labels=False): print("Det = {}".format(A.det())) D = matrix(eval(D)).rows() print("D:") print(D) G = point3d([v.list() for v in Dn(D,A,k)], size=8)#, opacity=.85) if labels: G += sum([text3d(str(v),v) for v in Dn(D,A,k)]) show(G, axes=False, frame=False)
Exploring Mandelbrot
Pablo Angulo
%cython
import numpy as np
cimport numpy as np
def mandelbrot_cython(float x0,float x1,float y0,float y1,
int N=200, int L=50, float R=3):
'''returns an array NxN to be plotted with matrix_plot
'''
cdef double complex c, z, I
cdef float deltax, deltay, R2 = R*R
cdef int h, j, k
cdef np.ndarray[np.uint16_t, ndim=2] m
m = np.zeros((N,N), dtype=np.uint16)
I = complex(0,1)
deltax = (x1-x0)/N
deltay = (y1-y0)/N
for j in range(N):
for k in range(N):
c = (x0+j*deltax)+ I*(y0+k*deltay)
z=0
h=0
while (h<L and
z.real**2 + z.imag**2 < R2):
z=z*z+c
h+=1
m[j,k]=h
return mimport pylab
x0_default = -2
y0_default = -1.5
side_default = 3.0
side = side_default
x0 = x0_default
y0 = y0_default
options = ['Reset','Upper Left', 'Upper Right', 'Stay', 'Lower Left', 'Lower Right']
@interact
def show_mandelbrot(option = selector(options, nrows = 2, width=8),
N = slider(100, 1000,100, 300),
L = slider(20, 300, 20, 60),
plot_size = slider(2,10,1,6),
auto_update = False):
global x0, y0, side
if option == 'Lower Right':
x0 += side/2
y0 += side/2
elif option == 'Upper Right':
y0 += side/2
elif option == 'Lower Left':
x0 += side/2
if option=='Reset':
side = side_default
x0 = x0_default
y0 = y0_default
elif option != 'Stay':
side = side/2
m=mandelbrot_cython(x0 ,x0 + side ,y0 ,y0 + side , N, L )
# p = (matrix_plot(m) +
# line2d([(N/2,0),(N/2,N)], color='red', zorder=2) +
# line2d([(0,N/2),(N,N/2)], color='red', zorder=2))
# time show(p, figsize = (plot_size, plot_size))
pylab.clf()
pylab.imshow(m, cmap = pylab.cm.gray)
pylab.savefig('mandelbrot.png')
Mandelbrot & Julia Interact with variable exponent
published notebook: https://cloud.sagemath.com/projects/19575ea0-317e-402b-be57-368d04c113db/files/pub/1201-1301/1299-Mandelbrot.sagews
Mandelbrot
by Harald Schilly
xxxxxxxxxxdef mandel_plot(expo = slider(-10,10,0.1,2), \ formula = ['mandel','ff'],\ iterations=slider(1,100,1,30), \ zoom_x = range_slider(-2,2,0.01,(-2,1)), \ zoom_y = range_slider(-2,2,0.01,(-1.5,1.5))): var('z c') f(z,c) = z^expo + c ff_m = fast_callable(f, vars=[z,c], domain=CDF) # messing around with fast_callable for i in range(int(iterations)/3): f(z,c) = f(z,c)^expo+c ff = fast_callable(f, vars=[z,c], domain=CDF) def mandel(z): c = z for i in range(iterations): z = ff_m(z,c) if abs(z) > 2: return z return z print('z <- z^%s + c' % expo) # calling ff three times, otherwise it fast_callable exceeds a recursion limit if formula is 'ff': func = lambda z: ff(ff(ff(z,z),z),z) elif formula is 'mandel': func = mandel complex_plot(func, zoom_x,zoom_y, plot_points=200, dpi=150).show(frame=True, aspect_ratio=1)
Julia
by Harald Schilly
xxxxxxxxxxdef julia_plot(expo = slider(-10,10,0.1,2), \ iterations=slider(1,100,1,30), \ c_real = slider(-2,2,0.01,0.5), \ c_imag = slider(-2,2,0.01,0.5), \ zoom_x = range_slider(-2,2,0.01,(-1.5,1.5)), \ zoom_y = range_slider(-2,2,0.01,(-1.5,1.5))): var('z') I = CDF.gen() f(z) = z^expo + c_real + c_imag*I ff_j = fast_callable(f, vars=[z], domain=CDF) def julia(z): for i in range(iterations): z = ff_j(z) if abs(z) > 2: return z return z print('z <- z^%s + (%s+%s*I)' % (expo, c_real, c_imag)) complex_plot(julia, zoom_x,zoom_y, plot_points=200, dpi=150).show(frame=True, aspect_ratio=1)
julia_plot(-7,30,0.5,0.5,(-1.5,1.5), (-1.5,1.5))
Sierpiński Triangle
by Eviatar Bach
xxxxxxxxxxdef sierpinski(N): '''Generates the Sierpinski triangle by taking the modulo-2 of each element in Pascal's triangle''' return [([0] * (N // 2 - a // 2)) + [binomial(a, b) % 2 for b in range(a + 1)] + ([0] * (N // 2 - a // 2)) for a in range(0, N, 2)]def _(N=slider([2 ** a for a in range(12)], label='Number of iterations', default=64), size=slider(1, 20, label='Size', step_size=1, default=9)): M = sierpinski(2 * N) matrix_plot(M, cmap='binary').show(figsize=[size, size])