Fractal pour faire et jouer avec Python

Fractal est un très bon matériel pédagogique pour les mathématiques et la programmation. Présentation de la figure fractale qui a été rayée après la pratique.

Je n'expliquerai pas la théorie.

Ensemble de Mandelbrot

Le code complet est ici.

En parlant de fractale, c'est un ensemble de Mandelbrot.

Définition $ z_0 = 0, z_ {n + 1} = z_n ^ 2 + point sur le plan complexe où la séquence complexe $ \ {z_n \} $ définie par C $ ne diverge pas en $ n \ vers \ infty $ L'ensemble de C $ est appelé l'ensemble de Mandelbrot.

Séquence de nombres complexes ici\\{z_n \\}Mais$n \to \infty Ne diverge pas|z_n|Maisn \to \infty Dit qu'il ne diverge pas. Aussi,|z_k| > 2$Rencontrerz_kMais存在する場合は\\{z_n \\}Maisいずれ発散することMais知られています。この命題の説明はIciC'est dedans.

Par conséquent, pour chaque $ C $, $ \ {z_n \} $ est calculé jusqu'à 100 éléments, et le jugement de convergence est effectué en fonction de l'existence ou non d'une valeur absolue supérieure à 2. Il existe différentes méthodes de coloration pour l'ensemble de Mandelbrot, mais cette fois, nous colorerons en fonction du nombre d'iters requis pour le jugement de convergence.

`mandelbrot_set.py`


solutions = np.array([2.7693, 0.11535-0.58974J, 0.11535+0.58974J])

def f(z, C):
    return z*z + C

def convergence_judgment(C):
    z = inital_z # 0.0
    i = 0
    while i < 100:
        z = f(z, C)   
        if abs(z) > 2:
            return i #Minimum de plus de 2|z_n|Renvoie n correspondant à.

        i += 1 
    return i #S'il ne diverge pas, renvoyez 100.

Dans la fonction principale, écrivez qu'un rectangle dont la longueur est epsilon centrée sur la coordonnée centrale center_position est subdivisé en N × N grilles, et le jugement de convergence est effectué dans chaque grille.

`mandelbrot_set.py`


N = 1000
inital_z = 0

epsilon = 1 #<2
center_position = [0, 0]

def main():
    grid = np.empty((N, N), dtype=np.uint8)
    points_x = np.linspace(center_position[0]-epsilon, center_position[0]+epsilon, N)
    points_y = np.linspace(center_position[1]-epsilon, center_position[1]+epsilon, N)

    for n, ReC in tqdm(enumerate(points_x)):
        for m, ImC in enumerate(points_y):
            C = ReC + ImC*1J
            grid[m, n] = convergence_judgment(C)

    show(grid)

Dans ce code, le temps de calcul était N = 1000 et il était d'environ 100 secondes. Cliquez sur l'image pour l'agrandir.


(0,0)	(-0.53 -0.613)	(-0.76,0.067)

De plus, bien que cela s'écarte de la définition de l'ensemble de Mandelbrot, si vous changez le premier terme $ z_0 $, vous pouvez revoir l'apparence caractéristique et c'est amusant! De plus, il est amusant de changer l'index 2 de z ^ 2 + C et de voir le modèle de z ^ n + C pour tout n!

Fractale de Newton

Le code complet est ici.

Par exemple, une équation cubique pour le nombre complexe $ z $

z^3-3z^2+z-1=0\tag{1}

A trois solutions. Comme je suis négligent, je peux trouver une solution en entrant ce qui suit dans Mathics.

Méthode de Newton avec n'importe quel point du plan complexe comme point de départ $ z_0 $

z_{n+1} = z_{n} - \frac{f(z_n)}{f'(z_n)}

Pensez à mettre à jour $ z $ avec pour trouver la solution de l'équation (1). Bien sûr, le côté gauche de l'équation (1) n'est pas monotone par rapport à $ z $, donc je ne sais pas si la méthode de Newton trouvera la solution. Le résultat de la méthode Newton pour un $ z_0 $ correctement sélectionné

Accédez à la solution 1
Accédez à la solution 2
Accédez à la solution 3
Je ne peux aller nulle part (il ne converge pas par la méthode Newton)

Vous obtiendrez l'un des résultats suivants. Pour chaque point de départ $ z_0 $ de la méthode de Newton, le point $ z_0 $ sur le plan complexe est coloré en fonction duquel le résultat de 1 à 4 ci-dessus est obtenu.

`newton_fractal.py`


def f(z):
    return z**3-3*z**2+z-1

def diff_f(z):
    return 3*z**2-6*z+1

def find_nearest_solution_idx_and_residual(z):
    nearest_solution_idx = abs(z-solutions).argmin()
    residual = abs(z - solutions[nearest_solution_idx])
    
    return nearest_solution_idx, residual 

def solve_f(inital_z):
    z = inital_z
    for n in range(max_step):
        z = z - f(z) / diff_f(z)
        nearest_solution_idx, residual = find_nearest_solution_idx_and_residual(z)
        if residual < residual_criteria: #convergence judgment
            return nearest_solution_idx + 2*n/max_step

    return len(solutions)

Ensuite, la fractale suivante apparaîtra. Le temps de calcul était N = 2500 et il était de 24 minutes.

C'est amusant de jouer avec différentes équations!

[1] [Système mécanique complexe et calculabilité Hiroyuki Inao](https://www.math.kyoto-u.ac.jp/insei/?plugin=attach&pcmd=open&file=Inou.pdf&refer=KINOSAKI%20SEMINAR%202009% 2Fattach)

Joint Shellpinsky

Le code complet est ici

Le joint du joint Shellpinsky est généré par l'algorithme suivant.

Déterminez correctement les coordonnées du triangle régulier (a1, a2, a3).
Reliez les milieux de chaque côté du triangle régulier (a1, a2, a3).
Sur les quatre triangles apparus à l'étape 2, faites 2 pour les trois triangles à l'exclusion du triangle central.

Le code source ci-dessous utilise des nombres complexes pour les coordonnées pour plus de commodité, mais ce n'est pas essentiel. triangle est une liste de trois éléments, chaque élément correspondant aux coordonnées complexes des sommets du triangle.

`sierpinski_gasket_fractals.py`


#Générez trois triangles à partir d'un triangle.
def make_3triangles_from_triangle(triangle):
    p0 = triangle[0]
    p1 = triangle[1]
    p2 = triangle[2]

    new_p2 = (m*p0 + n*p1) / (m+n)
    new_p0 = (m*p1 + n*p2) / (m+n)
    new_p1 = (m*p2 + n*p0) / (m+n)

    triangle0 = [p0, new_p1, new_p2]
    triangle1 = [new_p1, p2, new_p0]
    triangle2 = [new_p2, new_p0, p1]
    
    return triangle0, triangle1, triangle2 

def get_new_triangles(triangles):
    new_triangles = []
    for i, triangle in enumerate(triangles):
        triangle0, triangle1, triangle2 = make_3triangles_from_triangle(triangles[i])
        new_triangles.append(triangle0)
        new_triangles.append(triangle1)
        new_triangles.append(triangle2)

    return new_triangles

Puisque c'est un gros problème, trouvons la somme des circonférences et de la surface. Théoriquement, la circonférence devrait converger vers l'infini et l'aire vers zéro.

`sierpinski_gasket_fractals.py`


def get_circumference_and_area(triangles):
    distance_between_two_point = abs(triangles[-1][0] - triangles[-1][1])
    circumference = len(triangles) * distance_between_two_point
    area = math.sqrt(3)/4.0 * distance_between_two_point ** 2
    return circumference, area

Lorsque vous exécutez le programme, vous pouvez voir que la circonférence augmente et que la zone converge vers 0.

depth:0 num_triangles:3 circumference:1.50000 area:0.10825 time:0.000010 [sec]
depth:1 num_triangles:9 circumference:2.25000 area:0.02706 time:0.000118 [sec]
depth:2 num_triangles:27 circumference:3.37500 area:0.00677 time:0.000203 [sec]
depth:3 num_triangles:81 circumference:5.06250 area:0.00169 time:0.000349 [sec]
depth:4 num_triangles:243 circumference:7.59375 area:0.00042 time:0.000826 [sec]
depth:5 num_triangles:729 circumference:11.39062 area:0.00011 time:0.001953 [sec]
depth:6 num_triangles:2187 circumference:17.08594 area:0.00003 time:0.006996 [sec]
depth:7 num_triangles:6561 circumference:25.62891 area:0.00001 time:0.044693 [sec]
depth:8 num_triangles:19683 circumference:38.44336 area:0.00000 time:0.117012 [sec]
depth:9 num_triangles:59049 circumference:57.66504 area:0.00000 time:0.377997 [sec]
drawing...

J'ai fait une animation parce que c'était un gros problème.

J'ai changé un peu le code source et j'ai fabriqué un étrange joint Shelpinsky. Là encore, la circonférence diverge et la zone converge vers zéro.

Courbe de dragon

Le code complet est ici.

Comme le montre la figure ci-dessous, vous pouvez dessiner une image de courbe appelée courbe de dragon en ajoutant des images pivotées de 90 degrés autour du point final du segment de ligne à chaque étape.

Si le point $ (x, y) $ est tourné autour du point $ (C_x, C_y) $ de $ θ $, le point $ (x, y) $ sera

x_{new} = x\cos{\theta} - y\sin{\theta} + C_x - C_x\cos{\theta} + C_y\sin{\theta}  \\ 
y_{new} = x\sin{\theta} + y\cos{\theta} + C_y - C_x\sin{\theta} - C_y\cos{\theta}

Aller à. Rendez cela facile

x_{new} = f_x(x,y,C_{x},C_{y})  \\ 
y_{new} = f_y(x,y,C_{x},C_{y})

J'écrirai.

À une profondeur de $ j $, $ 2 ^ j $ de nouveaux points apparaîtront. Nous écrirons les $ 2 ^ j $ points nouvellement apparus sous la forme $ x_ {new} ^ {(k)} $, $ y_ {new} ^ {(k)} $. Cependant, $ k = 1, 2,…, 2 ^ j $. Puis, à la profondeur de $ j $, pour chaque $ k $

x_{new}^{(k)} = f_x(x_{2^j-k},y_{2^j-k},x_{2^j},y_{2^j})  \\ 
y_{new}^{(k)} = f_y(x_{2^j-k},y_{2^j-k},x_{2^j},y_{2^j})

Vous pouvez obtenir les coordonnées nouvellement ajoutées en calculant.

`DragonCurve.py`


def get_rotate_coordinate(x, y, Cx, Cy, theta=Pi/2):
    newx = x*math.cos(theta)-y*math.sin(theta)+Cx-Cx*math.cos(theta)+Cy*math.sin(theta)
    newy = x*math.sin(theta)+y*math.cos(theta)+Cy-Cx*math.sin(theta)-Cy*math.cos(theta)
    return newx, newy

def main():
    # Initial Condition
    xs = [1, 0]
    ys = [1, 1]

    for depth in range(max_depth):
        for k in range(2**depth):
            newx, newy = get_rotate_coordinate(xs[2**depth-k-1], ys[2**depth-k-1], Cx=xs[2**depth], Cy=ys[2**depth])
            xs.append(newx)
            ys.append(newy)
        print("depth:{} length:{}".format(depth, len(xs)-1))
    
    print("drawing..")
    draw(xs, ys)

Le chiffre ainsi obtenu est le suivant.


depth 5	depth 10	depth 20

Si vous décalez l'angle de rotation de $ \ theta = 90 ° $, la forme de la courbe du dragon sera continuellement transformée, ce qui est intéressant!