Fraktal zum Erstellen und Spielen mit Python

Fraktal ist ein sehr gutes Unterrichtsmaterial für Mathematik und Programmierung. Vorstellung der fraktalen Figur, die nach dem Training zerkratzt wurde.

Ich werde die Theorie nicht erklären.

Mandelbrot eingestellt

Der gesamte Code ist hier.

Apropos Fraktal, es ist ein Mandelbrot-Set.

Definition $ z_0 = 0, z_ {n + 1} = z_n ^ 2 + Punkt auf der komplexen Ebene, an dem die durch C $ definierte komplexe Folge $ \ {z_n \} $ bei $ n \ zu \ infty $ nicht divergiert Die Menge von C $ wird als Mandelbrot-Menge bezeichnet.

Komplexe Zahlenfolge hier\\{z_n \\}Aber$n \to \infty Geht nicht auseinander|z_n|Abern \to \infty Es heißt, dass es nicht divergiert. Ebenfalls,|z_k| > 2$Treffenz_kAber存在する場合は\\{z_n \\}Aberいずれ発散することAber知られています。この命題の説明はHierEs ist in.

Daher wird für jedes $ C $ $ \ {z_n \} $ für bis zu 100 Elemente berechnet, und die Konvergenzbeurteilung wird basierend darauf durchgeführt, ob es einen Absolutwert gibt, der 2 überschreitet oder nicht. Es gibt verschiedene Färbemethoden für das Mandelbrot-Set, aber dieses Mal färben wir entsprechend der Anzahl der für die Konvergenzbeurteilung erforderlichen Iter.

`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 #Mindestens mehr als 2|z_n|Geben Sie n entsprechend zurück.

        i += 1 
    return i #Wenn es nicht divergiert, geben Sie 100 zurück.

Schreiben Sie in der Hauptfunktion, dass ein Rechteck, dessen Länge epsilon zentriert auf der Mittelkoordinate center_position ist, in N × N Gitter unterteilt ist und in jedem Gitter eine Konvergenzbeurteilung durchgeführt wird.

`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)

In diesem Code betrug die Berechnungszeit N = 1000 und ungefähr 100 Sekunden. Klicken Sie auf das Bild, um es zu vergrößern.


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

Obwohl es von der Definition des Mandelbrot-Sets abweicht, können Sie beim Ändern des ersten Terms $ z_0 $ das charakteristische Erscheinungsbild wieder sehen und es macht Spaß! Es macht auch Spaß, den Index 2 von z ^ 2 + C zu ändern und das Muster von z ^ n + C für jedes n zu sehen!

Newton Fractal

Der gesamte Code lautet hier.

Zum Beispiel eine kubische Gleichung für die komplexe Zahl $ z $

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

Hat drei Lösungen. Ich bin fahrlässig. Wenn Sie also in Mathematik Folgendes eingeben, können Sie die Lösung finden.

Newton-Methode mit einem beliebigen Punkt auf der komplexen Ebene als Startpunkt $ z_0 $

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

Erwägen Sie, $ z $ mit zu aktualisieren, um die Lösung für Gleichung (1) zu finden. Natürlich ist die linke Seite von Gleichung (1) in Bezug auf $ z $ nicht eintönig, daher weiß ich nicht, ob die Newton-Methode die Lösung finden wird. Das Ergebnis der Newton-Methode für ein richtig ausgewähltes $ z_0 $

Fahren Sie mit Lösung 1 fort
Fahren Sie mit Lösung 2 fort
Fahren Sie mit Lösung 3 fort
Ich komme nicht weiter (es konvergiert nicht nach der Newton-Methode)

Sie erhalten eines der folgenden Ergebnisse. Für jeden Startpunkt $ z_0 $ der Newton-Methode wird der Punkt $ z_0 $ auf der komplexen Ebene gefärbt, gemäß dem das Ergebnis von 1 bis 4 oben erhalten wird.

`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)

Dann erscheint das folgende Fraktal. Die Berechnungszeit betrug N = 2500 und 24 Minuten.

Es macht Spaß, mit verschiedenen Gleichungen herumzuspielen!

[1] [Hiroyuki Inao Komplexes mechanisches System und Berechenbarkeit](https://www.math.kyoto-u.ac.jp/insei/?plugin=attach&pcmd=open&file=Inou.pdf&refer=KINOSAKI%20SEMINAR%202009% 2Fattach)

Shellpinsky Dichtung

Der gesamte Code ist hier

Die Dichtung für die Shellpinsky-Dichtung wird durch den folgenden Algorithmus erzeugt.

Bestimmen Sie die Koordinaten des regulären Dreiecks (a1, a2, a3) entsprechend.
Verbinden Sie die Mittelpunkte jeder Seite des regulären Dreiecks (a1, a2, a3).
Führen Sie von den vier in Schritt 2 angezeigten Dreiecken 2 für die drei Dreiecke mit Ausnahme des zentralen Dreiecks aus.

Der folgende Quellcode verwendet der Einfachheit halber komplexe Zahlen für die Koordinaten, dies ist jedoch nicht unbedingt erforderlich. Dreieck ist eine Liste mit drei Elementen, wobei jedes Element den komplexen Koordinaten der Eckpunkte des Dreiecks entspricht.

`sierpinski_gasket_fractals.py`


#Generieren Sie drei Dreiecke aus einem Dreieck.
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

Da es eine große Sache ist, lassen Sie uns die Summe der Umfänge und der Fläche ermitteln. Theoretisch sollte der Umfang gegen unendlich und die Fläche gegen null konvergieren.

`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

Wenn Sie das Programm ausführen, können Sie sehen, dass der Umfang zunimmt und die Fläche gegen 0 konvergiert.

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...

Ich habe eine Animation gemacht, weil es eine große Sache war.

Ich habe den Quellcode ein wenig geändert und eine seltsame Shelpinsky-Dichtung hergestellt. Wiederum divergiert der Umfang und die Fläche konvergiert gegen Null.

Drachenkurve

Der gesamte Code lautet hier.

Wie in der folgenden Abbildung gezeigt, können Sie ein Kurvenbild zeichnen, das als Drachenkurve bezeichnet wird, indem Sie bei jedem Schritt Bilder hinzufügen, die um 90 Grad um den Endpunkt des Liniensegments gedreht sind.

Wenn der Punkt $ (x, y) $ um $ θ $ um den Punkt $ (C_x, C_y) $ gedreht wird, ist der Punkt $ (x, y) $

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}

Gehe zu. Mach das einfach

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

Ich werde schreiben.

In einer Tiefe von $ j $ erscheinen $ 2 ^ j $ neue Punkte. Wir werden die neu erscheinenden $ 2 ^ j $ -Punkte als $ x_ {new} ^ {(k)} $, $ y_ {new} ^ {(k)} $ schreiben. $ K = 1, 2,…, 2 ^ j $. Dann in der Tiefe von $ j $ für jedes $ 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})

Sie können die neu hinzugefügten Koordinaten durch Berechnung erhalten.

`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)

Die auf diese Weise erhaltene Figur ist wie folgt.


depth 5	depth 10	depth 20

Wenn Sie den Drehwinkel von $ \ theta = 90 ° $ verschieben, wird die Form der Drachenkurve kontinuierlich geändert, was interessant ist!