Schnelles rechnen, Matrix, Konvergenz, Julia-Fraktale

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
Ctrl-Z
User
Beiträge: 12
Registriert: Freitag 1. Mai 2015, 03:55

Hi liebe Leute,

ich möchte gerade Julia-Fraktale ausrechnen. Das Problem ist folgendes: Ich rechne die Iterationen für ein vorgebenes Gitter, zum Beispiel 500 * 500. Nun gibt es natürlich das Problem, dass viele Gitterpunkte divergieren, die schießen schnell ins unendliche.
Man muss also definitiv eine Kontrollfunktion einbauen, die einem die Begrenztheit ermöglicht. Nun habe ich mir was überlegt (das Bild ist auch sehr schön anzuschauen!), und zwar benutze ich Matrix aus Nullen und Einsen, die angibt, ob das Grenzkriterium erreicht ist. Mit der multipliziere ich dann, und dank der genialen Struktur von Python ergibt sich ein Ergebnis, was niemals divergiert, juhu!
Blöde nur: Es ist mathemathisch nicht korrekt, da genullte Positionen im Gitter ständig neu iteriert werden! Dadurch wird das Bild hübscher und bunter, aber ist nicht korrekt und ist natürlich viel langsamer, da in jeder Iteration jeder Punkt neu berechnet wird. Hier mein Quellcode:

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt

def main():
    # complex constant
    cr = -0.4
    ci = +1.0
    # real domain
    r0=-2.0
    r1=2.0
    # complex domain
    i0=-2.0
    i1=2.0
    # resolution
    divisions = 500
    # cutoff
    limit = 2
    #number of iterations
    niter = 30

    # define mesh
    rr=np.linspace(r0,r1,divisions)
    ii=np.linspace(i0,i1,divisions)
    rv,iv = np.meshgrid(rr,ii)
    #print rv,iv

    #calculation
    for k in range(niter):
        rv1 = rv*rv - iv*iv + cr # real part
        iv1 = 2 * iv * rv + ci # imaginary part
        rvbol = (np.abs(rv1) <= limit) # boolean: 1 if True
        ivbol = (np.abs(iv1) <= limit)
        rv = rv1 *rvbol # set to zero if above limit
        iv = iv1 *ivbol
        del rvbol, ivbol, rv1, iv1


    # contour levels
    levels = np.linspace(0,2,200)
    # plot graph
    zabs = np.abs(np.sqrt(rv**2+iv**2))
    plt.contourf(rr,ii,zabs,levels=levels)
    plt.show()
    pass

if __name__ == '__main__':
    main()

Der Algorhitmus soll aber nur Punkte berücksichtigen, die nicht durch Divergenz eh schon gebrandmarkt sind. Dabei möchte ich die numerische Stärke von Numpy ausnutzen, dass erscheint mir eigentlich sinnvoll. Ich denke gerade über Sparse Matrizen nach, bringt das was?
Das wäre doch super, wenn ich noch Punkte multipliziere, die relevant sind!
Was meint ihr?

LG,
Jens
Ctrl-Z
User
Beiträge: 12
Registriert: Freitag 1. Mai 2015, 03:55

Das Thema ist glaube ich durch, ich benutze eine Boolean-Matrix zur Registrierung von Grenzüberschreitungen und eine klassische Doppelschleife über die Koordinaten nebst einer über die Iterationen, wunderbar! Ganz einfach! 8)

Ist mathemaisch korrekt und performanter, da sofort geprüft wird, ob die Grenze schon überschritten wurde. Blöd nur, dass da Numpy zu kurz kommt... vielleicht geht das ja noch zig mal schneller? Danke für Ideen, aber ich bin erstmal glücklich.

Code: Alles auswählen

calc = 0

divergedij = np.zeros([yres,xres],dtype=bool)
z = np.zeros([yres,xres],dtype=float)


for itr in range(iterations):
    for i in range(xres):
        for j in range(yres):
            if not divergedij[j,i]:
                calc += 1
                xm1 = xm[j,i]**2 - ym[j,i]**2 + rconst
                ym1 = 2 * xm[j,i] * ym[j,i] + iconst
                zabs = xm1**2 + ym1**2
                if (zabs) > limit:
                    divergedij[j,i] = True
                    z[j,i]=-1
                else:
                    xm[j,i] = xm1
                    ym[j,i] = ym1
                    z[j,i] = zabs


LG,
Jens
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

Die schnellste Mandelbrot Berechnung mit Numpy die ich kenne:
https://thesamovar.wordpress.com/2009/0 ... and-numpy/
a fool with a tool is still a fool, www.magben.de, YouTube
Antworten