Hintergrund

Seit dem letzten Post sind ungefähr zwei Jahre vergangen. Wir haben im vorherigen Artikel verschiedene Sprachen implementiert und die Berechnungsgeschwindigkeiten verglichen. Die Motivation ist, dass es eine versteckte Linie gibt, die die Berechnungsgeschwindigkeit durch willkürliche Genauigkeitsberechnung vergleichen möchte. Für die Berechnung der willkürlichen Genauigkeit ist die Implementierung von GNU MP (GMP), MPFR (sowohl C als auch C ++) von großer Bedeutung. ， Die Wrapper und Erweiterungen dieser Pakete sind in einigen Sprachen verfügbar. In Python ist es in mpmath verfügbar (ich kenne Decimal der Standardbibliothek nicht), und Julia und Mathematica (Wolfram-Sprache) unterstützen standardmäßig Operationen mit mehreren Längen. Beide scheinen GMP und MPFR zu verwenden, daher ist das ultimative Ziel zu wissen, wie viel Aufwand und einfache Implementierung der Codierung sind (diesmal nicht erreicht).

Übrigens wird die 4x-Genauigkeitsberechnung ab gcc 4.6 unterstützt, und die 4x-Genauigkeitsberechnung sollte in C, C ++ und Fortran möglich sein, wird hier jedoch nicht behandelt. (In C benötigte der Gleitkommatyp mit vierfacher Genauigkeit eine Deklaration wie __float128, aber was ist jetzt ...)

Beispiele, die eine beliebige Genauigkeit erfordern

(Eindimensional) [Bäckerkarte](https://ja.wikipedia.org/wiki/%E3%83%91%E3%82%A4%E3%81%93%E3%81%AD % E5% A4% 89% E6% 8F% 9B) $ B(x) = \begin{cases} 2x,& 0 \le x \le 1/2 \\\\ 2x-1,& 1/2 < x \le 1 \end{cases} $ Betrachten Sie $ x \ in [0,1] $ als Anfangsbedingung und die Orbitalsequenz $ \ {x, B (x), B ^ {(2)} = B \ circ, die durch iterative Synthese von $ B $ erhalten wird In Anbetracht von B (x), \ ldots, B ^ {(n)} (x) \} $ verhält sich dies chaotisch (das Brot ist durch die Operation von B gleichmäßig verbunden).

Schreiben Sie vorerst (ungefähr) Code, der Iterationen mit geeigneten Anfangsbedingungen bis in den dritten Stock zeichnet.

import numpy as np
import matplotlib.pyplot as plt
import time

def B(x):
    return 2*x if x <= 1/2 else (2*x - 1)

fig, ax = plt.subplots(1, 1, figsize=(5,5))
ax.set_xlabel("x")
ax.set_ylabel("B(x)")

x = np.linspace(0, 1, 300)
ax.set_ylim(-0.02,1.02)
ax.set_xlim(-0.02,1.02)
ax.plot(x,2*x,"-y")
ax.plot(x,2*x-1,"-y")
ax.plot(x,x,"--k")

x = (np.sqrt(5) - 1)/2 

web = [np.array([x]),
        np.array([0])]
traj = np.array([x])

nmax = 3

bboxprop = dict(facecolor='w', alpha=0.7)
ax.text(x+0.05,0, r"$x=%f$" % x, bbox=bboxprop)
a = time.time()
for n in range(1,nmax+1):
    y = B(x)
    web[0] = np.append(web[0], x)
    web[1] = np.append(web[1], y)
    traj = np.append(traj, y)
    ax.text(x+0.02, y, r"$f^{(%d)}(x)$" % n, bbox=bboxprop)

    x = y
    if n != nmax:
        ax.text(y+0.02, y, r"$x=f^{(%d)}(x)$" % n, bbox=bboxprop)
        
        web[0] = np.append(web[0], x)
        web[1] = np.append(web[1], y)        

print(time.time() - a, "[sec.]")
ax.plot(web[0][:], web[1][:], "-o")
ax.plot(web[0][0], web[1][0], "or")
ax.plot(web[0][-1], web[1][-1], "or")
ax.grid()
print(traj)
plt.savefig("web.png ")
plt.show()

Um die Orbitalsequenz von $ \ {x, B (x), B \ circ B (x), B \ circ B \ circ B (x) \} $ zu erhalten

$ python baker.py 
0.0013127326965332031 [sec.]
[0.61803399 0.23606798 0.47213595 0.94427191]

Ich denke auch, dass die Abbildung (Webdiagramm oder Phasenporträt) die Bedeutung der iterativen Synthese visuell verstehen kann (Startpunkt und Endpunkt sind rot eingekreist).

Was ist nun, wenn wir die Anzahl der Iterationen auf 50 setzen?

nmax = 3

nmax = 50

Und Renn

$ python baker-2.py
0.0012233257293701172 [sec.]
[0.61803399 0.23606798 0.47213595 0.94427191 0.88854382 0.77708764
 0.55417528 0.10835056 0.21670112 0.43340224 0.86680448 0.73360896
 0.46721792 0.93443584 0.86887168 0.73774336 0.47548671 0.95097343
 0.90194685 0.8038937  0.60778741 0.21557482 0.43114964 0.86229928
 0.72459856 0.44919711 0.89839423 0.79678845 0.59357691 0.18715382
 0.37430763 0.74861526 0.49723053 0.99446106 0.98892212 0.97784424
 0.95568848 0.91137695 0.82275391 0.64550781 0.29101562 0.58203125
 0.1640625  0.328125   0.65625    0.3125     0.625      0.25
 0.5        1.         1.        ]

Das Ergebnis von wurde erhalten. Da $ B (1) per Definition = 1 $ ist, wird eine ausreichend lange iterative Synthese $ B ^ {(n)} (x) $ mit der Anfangsbedingung $ x = (\ sqrt {5} -1) / 2 $ auf 1 gesetzt. Wird voraussichtlich von der numerischen Berechnung konvergieren (Pan mischt sich nicht gleichmäßig ...)

Diese auf numerischer Berechnung basierende Vorhersage wurde jedoch aufgrund unzureichender Berechnungsgenauigkeit fälschlicherweise abgeleitet. Dieses System verhält sich für fast alle Anfangspunkte chaotisch, mit Ausnahme spezieller Anfangsbedingungen ($ x = 0,1 / 2,1 $ usw.), wodurch $ x \ in [0,1] $ dicht wird. ergänze.

import numpy as np
import matplotlib.pyplot as plt
import mpmath
import time
mpmath.mp.dps = 100

def B(x):
    return 2*x if x <= 1/2 else (2*x - 1)

fig, ax = plt.subplots(1, 1, figsize=(5,5))
ax.set_xlabel("x")
ax.set_ylabel("B(x)")

x = np.linspace(0, 1, 300)
ax.set_ylim(-0.02,1.02)
ax.set_xlim(-0.02,1.02)
ax.plot(x,2*x,"-y")
ax.plot(x,2*x-1,"-y")
ax.plot(x,x,"--k")

x = (np.sqrt(5) - 1)/2 

web = [np.array([x]),
        np.array([0])]
traj = np.array([x])

x = (mpmath.sqrt(5) - 1)/2 
y = B(x)

web = [np.array([x]),
        np.array([0])]

traj = np.array([x])

nmax = 50
a = time.time()
for n in range(1,nmax+1):
    y = B(x)
    web[0] = np.append(web[0], x) 
    web[1] = np.append(web[1], y)
    traj = np.append(traj, y)    
    x = y
    if n!=nmax:
        web[0] = np.append(web[0], x)
        web[1] = np.append(web[1], y)    
print(time.time() - a, "[sec.]")

ax.plot(web[0], web[1], "-o")
ax.plot(web[0][0], web[1][0], "or")
ax.plot(web[0][-1], web[1][-1], "or")
ax.grid()

mpmath.nprint(traj.tolist())
plt.savefig("web-mp1.png ")
plt.show()

Die wesentliche Änderung ist

import mpmath
mpmath.mp.dps = 100

Und Anfangsbedingungen

x = (mpmath.sqrt(5) - 1)/2

Nur. mpmath.mp.dps = 100 gibt an, dass die Auswertung im formalen Teil mit einer Genauigkeit von 100 Stellen durchgeführt werden soll. (Standardmäßig ist die doppelte Genauigkeit (15) angegeben.) Wenn print (traj) verwendet wird, werden unverändert 100 Ziffern ausgegeben. Konvertieren Sie daher den Traj von mpmath einmal in eine Liste und verwenden Sie nprint, um die Anzahl der angezeigten Ziffern zu verringern. Das Ausgabeergebnis ist

$ python baker_mp.py 
0.0019168853759765625 [sec.]
[0.618034, 0.236068, 0.472136, 0.944272, 0.888544, 0.777088, 0.554175, 0.108351, 0.216701, 0.433402, 0.866804, 0.733609, 0.467218, 0.934436, 0.868872, 0.737743, 0.475487, 0.950973, 0.901947, 0.803894, 0.607787, 0.215575, 0.43115, 0.862299, 0.724599, 0.449197, 0.898394, 0.796788, 0.593577, 0.187154, 0.374308, 0.748615, 0.49723, 0.994461, 0.988921, 0.977842, 0.955685, 0.911369, 0.822739, 0.645478, 0.290956, 0.581912, 0.163824, 0.327647, 0.655294, 0.310589, 0.621177, 0.242355, 0.48471, 0.96942, 0.93884]

Im Gegensatz zum Berechnungsergebnis mit doppelter Genauigkeit konvergiert es nicht in 50 Schritten gegen 1. Die Berechnungsgeschwindigkeit ist etwas langsamer, liegt jedoch im zulässigen Bereich. Lassen Sie uns die Zeit länger entwickeln Lassen Sie uns die Zeit mit nmax = 350 entwickeln Die Ausgabe wird auch entsprechend der Anzahl der Schritte länger, wenn Sie also die letzten 30 Schritte oder höher anzeigen

mpmath.nprint(traj.tolist()[-30:])

$ python baker_mp.py 
0.013526201248168945 [sec.]
[0.0256348, 0.0512695, 0.102539, 0.205078, 0.410156, 0.820313, 0.640625, 0.28125, 0.5625, 0.125, 0.25, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

Es wird. Selbst bei der Berechnung von 100 Stellen des falschen Teils scheint sich die Anzahl der Schritte 1 zu nähern, wenn sie lang genug ist. Was ist, wenn wir die Genauigkeit weiter verbessern?

mpmath.mp.dps = 200

Ausführen als.

$ python baker_mp.py 
0.013622760772705078 [sec.]
[0.0256048, 0.0512095, 0.102419, 0.204838, 0.409676, 0.819353, 0.638705, 0.27741, 0.55482, 0.109641, 0.219282, 0.438564, 0.877127, 0.754255, 0.50851, 0.017019, 0.0340381, 0.0680761, 0.136152, 0.272304, 0.544609, 0.0892179, 0.178436, 0.356872, 0.713743, 0.427487, 0.854974, 0.709947, 0.419894, 0.839788]

Diesmal schien es sich nicht allmählich 1 zu nähern. Um eine langfristige Zeitentwicklung genau zu erhalten, ist es notwendig, die Genauigkeit zu erhöhen.

Vergleich der Rechengeschwindigkeit durch mehrere Implementierungen

Da die Implementierung des obigen Beispiels np.append verwendet, wird die Berechnungsgeschwindigkeit extrem langsam. In meiner Umgebung

mpmath.mp.dps = 100000
nmax = 10000

Bei Ausführung mit

$ python baker_mp.py 
1.895829439163208 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]

Es dauert etwas weniger als 2 Sekunden. Obwohl das Anhängen praktisch ist, ist es langsam, da es sich um eine Operation handelt, bei der die erneute Erfassung eines Arrays wiederholt wird (Wiederholen von Malloc und Kopieren). Bereiten Sie die erforderlichen Sequenzen im Voraus vor und schließen Sie sie nur mit der Zuweisungsoperation ab. Es sollte schneller sein.

Da die Klasse der Matrix in mpmath definiert ist, verwenden wir sie. Da es sich jedoch um eine Python-Liste handelt, die einem eindimensionalen Array entspricht, vergleichen Sie den Fall der Liste mit dem Fall der numpy.array. Der folgende Code zeigt nur die geänderten Teile


nmax = 10000
web = mpmath.zeros(2*nmax+1,2)
web[0,0] = x
web[0,1] = 0
traj = [mpmath.mpf("0")]*(nmax+1)
traj[0] = x

a = time.time()

for n in range(1,nmax+1):
    y = B(x)
    web[2*n - 1, 0] = x
    web[2*n - 1, 1] = y
    traj[n] = y
    x = y            
    if n!=nmax:
        web[2*n, 0] = x
        web[2*n, 1] = y

print(time.time()-a, "[sec.]")
mpmath.nprint(traj[-30:])

$ python baker_mp1.py 
0.355421781539917 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]

Versuchen Sie als nächstes, traj mit numpy.array abzudecken

traj = np.array([mpmath.mpf("0")]*(nmax+1))

dtype wird zum Objekt. Alle Ziffern werden ausgegeben, es sei denn, die Ausgabe wird einmal an die Liste zurückgegeben.

mpmath.nprint(traj[-30:].tolist())

Das Ausführungsergebnis scheint fast das gleiche zu sein.

$ python baker_mp2.py 
0.36023831367492676 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]

Als nächstes ist es kompliziert, aber füge numpy.array in die Liste ein.

web = [np.array([mpmath.mpf(0)]*2*(nmax+1)),
        np.array([mpmath.mpf(0)]*2*(nmax+1))]
traj = np.array([mpmath.mpf(0)]*(nmax+1))
traj[0] = x
        
a = time.time()

for n in range(1,nmax+1):
    y = B(x)
    web[0][2*n - 1] = x
    web[1][2*n - 1] = y
    traj[n]  = y
    x = y
    if n != nmax:
        web[0][2*n] = x
        web[1][2*n] = y    

print(time.time()-a, "[sec.]")
mpmath.nprint(traj[-30:].tolist())

Das Ausführungsergebnis ist etwas schneller.

$ python baker_mp3.py 
0.2923922538757324 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]

Mit numpy.array

web = np.zeros((2*nmax+1,2),dtype="object")
web[0,0] = x
web[0,1] = 0
traj = np.zeros(nmax+1, dtype="object")
traj[0] = x

a = time.time()
for n in range(1,nmax+1):
    y = B(x)
    web[2*n - 1, 0] = x
    web[2*n - 1, 1] = y
    traj[n] = y
    x = y            
    if n!=nmax:
        web[2*n, 0] = x
        web[2*n, 1] = y


print(time.time()-a, "[sec.]")
mpmath.nprint(traj[-30:].tolist())

Auch wenn sich nicht viel ändert.

$ python baker_mp4.py 
0.2934293746948242 [sec.]
[0.170384, 0.340767, 0.681535, 0.36307, 0.72614, 0.452279, 0.904559, 0.809118, 0.618235, 0.236471, 0.472941, 0.945882, 0.891765, 0.78353, 0.56706, 0.13412, 0.26824, 0.536479, 0.0729585, 0.145917, 0.291834, 0.583668, 0.167335, 0.334671, 0.669341, 0.338683, 0.677365, 0.35473, 0.70946, 0.418921]

Speichern Sie die Daten mit doppelter Genauigkeit.

web = np.zeros((2*nmax+1,2),dtype=np.float)
web[0,0] = float(x)
web[0,1] = 0
traj = np.zeros(nmax+1, dtype=np.float)
traj[0] = float(x)

a = time.time()
for n in range(1,nmax+1):
    y = B(x)
    web[2*n - 1, 0] = float(x)
    web[2*n - 1, 1] = float(y)
    traj[n] = float(y)
    x = y            
    if n!=nmax:
        web[2*n, 0] = float(x)
        web[2*n, 1] = float(y)

print(time.time()-a, "[sec.]")
print(traj[-30:])

Die Geschwindigkeit ändert sich fast nicht.

$ python baker_mp5.py 
0.3071897029876709 [sec.]
[0.17038373 0.34076746 0.68153493 0.36306985 0.72613971 0.45227941
 0.90455883 0.80911766 0.61823531 0.23647062 0.47294124 0.94588249
 0.89176498 0.78352995 0.5670599  0.13411981 0.26823961 0.53647923
 0.07295846 0.14591691 0.29183383 0.58366766 0.16733532 0.33467064
 0.66934128 0.33868256 0.67736512 0.35473024 0.70946048 0.41892095]

Fazit

Ich habe verschiedene Implementierungsmethoden ausprobiert, bin jedoch zu dem Schluss gekommen, dass sich die Geschwindigkeit von Python grundsätzlich nicht ändert, wenn nicht das Anhängen verwendet wird.

Dieses Mal wurde nur eine Anfangsbedingung angestrebt, aber bei der Berechnung der Statistik ist eine Stichprobe der Anfangsbedingung erforderlich. Ich sollte es auch in diesem Fall tun, aber dieses Mal werde ich mein Bestes geben.

Referenz

mpmath
Julia:BigFloats and BigInts
Über die interne Implementierung von Mathematica (numerische Werte und verwandte Funktionen)
Einführung in Chaos Mechanical Systems 2. Auflage Originalautor: Robert L. Devaney (Übersetzer: Kenichi Goto, Neuauflage Übersetzer: Koji Kokubun, Yutaka Ishii, Shunsaku Arai, Masashi Kisaka) Kyoritsu Publishing

Vergleich der Berechnungsgeschwindigkeit durch Implementierung von Python mpmath (willkürliche Genauigkeitsberechnung) (Hinweis)

Hintergrund

Beispiele, die eine beliebige Genauigkeit erfordern

Vergleich der Rechengeschwindigkeit durch mehrere Implementierungen

Fazit

Referenz