Differences between revisions 23 and 24
Revision 23 as of 2010-11-22 05:28:04
Size: 8161
Editor: Eviatar
Comment:
Revision 24 as of 2010-11-22 05:32:45
Size: 8235
Editor: Eviatar
Comment:
Deletions are marked like this. Additions are marked like this.
Line 281: Line 281:
    '''Generates Sierpinski's Triangle to N iterations using the Rule 90 Elementary Cellular Automaton'''     '''Generates the Sierpinski's Triangle fractal to N iterations using the Rule 90 elementary cellular automaton. N is in powers of 2 because these produce "whole" triangles.'''

Sage Interactions - Fractal

goto interact main page

Mandelbrot's Fractal Binomial Distribution

def 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 v

html("<h1>Mandelbrot's Fractal Binomial Measure</h1>")

@interact
def _(mu0=(0.3,(0.0001,0.999)), k=(3,(1..14)), thickness=(1.0,(0.1,0.2,..,1.0))):
    v = muk_plot(mu0,k)
    line(v,thickness=thickness).show(xmin=0.5, xmax=0.5, ymin=0, figsize=[8,3])

binomial.png

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])?)

A = matrix([[1,1],[-1,1]])
D = [vector([0,0]), vector([1,0])]

@interact
def f(A = matrix([[1,1],[-1,1]]), D = '[[0,0],[1,0]]', k=(3..17)):
    print "Det = ", 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)

1.png

Demonstrating that the Twin Dragon Matrix is likely to yield a Tiling of a Compact Interval of R^2 as k->infinity (It does!)

A = matrix([[1,1],[-1,1]])
D = [vector([0,0]), vector([1,0])]

@interact
def f(A = matrix([[1,1],[-1,1]]), D = '[[0,0],[1,0]]', k=(3..17)):
    print "Det = ", 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)

2.png

Now in 3D

A = 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 ans
    
@interact
def f(A = matrix([[0,0,2],[1,0,1],[0,1,-1]]), D = '[[0,0,0],[1,0,0]]', k=(3..15), labels=True):
    print "Det = ", 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)

3.png

4.png


CategoryCategory

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 m

import 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':
        pass
    else:
        side = side/2
    
    time 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)
    time pylab.savefig('mandelbrot.png')

mandelbrot_cython.png

Mandelbrot & Julia Interact with variable exponent

published notebook: http://sagenb.org/pub/1299/

Mandelbrot

by Harald Schilly

@interact
def mandel_plot(expo = slider(-10,10,0.1,2), \
      formula = list(['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)

mandel-interact-02.png

Julia

by Harald Schilly

@interact
def 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-interact-01.png

julia_plot(-7,30,0.5,0.5,(-1.5,1.5), (-1.5,1.5))

julia-fractal-exponent--7.png

Sierpinski's Triangle

by Eviatar Bach

%python

from numpy import zeros

def sierpinski(N):
    '''Generates the Sierpinski's Triangle fractal to N iterations using the Rule 90 elementary cellular automaton. N is in powers of 2 because these produce "whole" triangles.'''
    M=zeros( (N,2*N+1), dtype=int)
    M[0,N]=1
    rule=[0, 1, 0, 1, 1, 0, 1, 0]
    
    for j in range(1,N):
        for k in range(N-j,N+j+1):
            l = 4*M[j-1,k-1] + 2*M[j-1,k] + M[j-1,k+1]
            M[j,k]=rule[ l ]
    return M
    
@interact
def _(N=slider([2**a for a in range(0, 12)], label='Number of iterations',default=128), size = slider(1, 20, label= 'Size', step_size=1, default=9 )):
    M = sierpinski(N)
    plot_M = matrix_plot(M, cmap='binary')
    plot_M.show( figsize=[size,size])

sierpinski.png

interact/fractal (last edited 2019-04-06 16:11:28 by chapoton)