[PYTHON] Beschleunigung der numerischen Berechnung mit NumPy / SciPy: Anwendung 2

Ziel

Es ist für Anfänger von Python und NumPy geschrieben. Es kann besonders nützlich sein für diejenigen, die eine numerische Berechnung physikalischer Systeme anstreben, die "Ich kann C-Sprache verwenden, habe aber kürzlich Python gestartet" oder "Ich verstehe nicht, wie man wie Python schreibt" entspricht. Vielleicht.

Außerdem kann es aufgrund meines mangelnden Studiums zu falschen Beschreibungen kommen. Bitte verzeihen Sie mir.

Zusammenfassung

Beschleunigung der numerischen Berechnung mit NumPy: Basics Beschleunigung der numerischen Berechnung mit NumPy / SciPy: Anwendung 1

Dies ist eine Fortsetzung von. Ich werde der grundlegenden numerischen Berechnungsmethode folgen, aber diesmal enthält sie auch ein wenig fortgeschrittenen Inhalt.

Algebra / Transzendentale Gleichung

Algebraische Gleichungen sind sogenannte handlösbare Gleichungen. Transzendentale Gleichungen sind ziemlich übertriebene Namen, aber sie sind Gleichungen, die mit algebraischen Methoden nicht gelöst werden können.

\sin(x) = \frac{x}{2}

Dies ist ein Kind. Diese Gleichung lautet: "Sie benötigen $ \ frac {x} {2} $, um $ \ sin (x) $ zu kennen, und $, um $ \ frac {x} {2} $ zu kennen. Es wird auch als "selbstkonsistente Gleichung" bezeichnet, da es die Struktur "\ sin (x) $ ist erforderlich" hat. Natürlich kann es numerisch gelöst werden, auch wenn es nicht algebraisch gelöst werden kann. Die Algorithmen umfassen die "Newton-Methode" und die "Dichotomie". Ich werde die Details des Algorithmus weglassen, aber es ist nicht so schwierig. Wenn Sie es nicht wissen, probieren Sie es aus. Es ist unvermeidlich, eine for-Schleife zu verwenden, da dies ähnlich ist. Die Newton-Methode konvergiert grundsätzlich schnell, aber eine Funktion mit guten Eigenschaften ist erforderlich, sofern es sich um eine differenzierbare Funktion handelt. Andererseits kann die Dichotomie immer gefunden werden, wenn es eine Lösung im angegebenen Intervall gibt, aber die Anzahl der Iterationen bis zur Konvergenz groß ist.

Wie Sie vielleicht bereits bemerkt haben, ** sind diese bereits in einem Modul namens scipy.optimize implementiert. ** For-Schleifen sind in der Funktion vollständig und laufen schnell. ..

Implementierung

Finden Sie heraus, wo die Lösung für die obige Gleichung ungefähr liegt:

figure5.png

Sie können sehen, dass es sich um $ x \ in (1.5, 2.0) $ handelt. Ich werde dies in einer Dichotomie lösen:

from scipy.optimize import bisect
def f(x):
	return np.sin(x) - x/2
	
bisect(f, 1.5, 2) 
>>> 1.895494267034337

Es ist sehr einfach. Wenn das Intervall zwei oder mehr Lösungen enthält, wissen Sie nicht, welche konvergieren wird. Stellen Sie daher sicher, dass Sie eine haben. Natürlich können Sie auch algebraische Gleichungen verwenden.

Fourier-Umwandlung

Apropos Fourier-Umwandlung, wir sind mit Signalverarbeitung usw. vertraut, aber ich werde auch dieses Mal mit der Geschichte der Physik fortfahren. Bitte haben Sie Verständnis dafür, dass es etwas mehr mathematische Gespräche geben wird. Erinnern Sie sich an die Diffusionsgleichung, die im vorherigen Artikel behandelt wurde. Bitte gib mir:

\frac{\partial}{\partial t}f(x, t) = \frac{\partial^2}{\partial x^2}f(x, t)

Dieses Kind könnte mit der Fourier-Transformation gelöst werden. Mal sehen, wie. Wir definieren die Fourier-Transformation und die inverse Fourier-Transformation wie folgt:

\tilde{f}(k, t) = \frac{1}{\sqrt{2\pi}}\int dx e^{-ikx}f(k, t) = {\cal F}[f(x, t)]\\
f(x, t) = \frac{1}{\sqrt{2\pi}}\int dk e^{ikx}\tilde{f}(k, t) = {\cal F}^{-1}[\tilde{f}(k, t)]

Lassen Sie uns diese Definition in die Diffusionsgleichung einsetzen, die es uns ermöglicht, eine $ x $ -Differenzierung durchzuführen:

\frac{\partial}{\partial t}\left(\frac{1}{\sqrt{2\pi}}\int dk e^{ikx}\tilde{f}(k, t)\right) = \frac{\partial^2}{\partial x^2}\left(\frac{1}{\sqrt{2\pi}}\int dk e^{ikx}\tilde{f}(k, t)\right)\\
\frac{\partial}{\partial t}\left(\frac{1}{\sqrt{2\pi}}\int dk e^{ikx}\tilde{f}(k, t)\right) = -\left(\frac{1}{\sqrt{2\pi}}\int dk k^2e^{ikx}\tilde{f}(k, t)\right)\\

Wenn Sie die $ k $ -Integration entfernen, wird sie außerdem zu einer Differentialgleichung erster Ordnung für $ t $, sodass sie leicht gelöst werden kann:

\frac{\partial}{\partial t}\tilde{f}(k, t) = -k^2\tilde{f}(k, t)\\
\tilde{f}(k, t) = e^{-k^2 t}\tilde{f}(k, 0)

Und ich werde eine inverse Fourier-Umwandlung davon machen:

\frac{1}{\sqrt{2\pi}}\int dk e^{ikx}\tilde{f}(k, t) = \frac{1}{\sqrt{2\pi}}\int dk e^{ikx}e^{-k^2 t}\tilde{f}(k, 0)\\
f(x, t) = \frac{1}{\sqrt{2\pi}}\int dk e^{ikx}e^{-k^2 t}\tilde{f}(k, 0)

Ich bin nicht sicher, ob es bei $ \ tilde {f} $ bleibt, also werde ich es auf $ f $ korrigieren:

f(x, t) = \frac{1}{2\pi}\int dk e^{ikx}e^{-k^2 t}\int dx' e^{-ikx'}f(x', 0)

Wenn die Gleichung für die Anfangsbedingung $ f (x, 0) $ bekannt ist und analytisch integriert werden kann, kann die analytische Lösung [^ 1] erhalten werden. Wenn sie nicht analytisch integriert werden kann. Immerhin ist es nicht numerisch gelöst.

Eine einfache Beschreibung dieses Flusses lautet nun wie folgt:

f(x, t) = {\cal F}^{-1}[e^{-k^2 t}{\cal F}[f(x, 0)]]

Nachdem Sie $ f (x, 0) $ in Fourier konvertiert haben, multiplizieren Sie es mit $ e ^ {-k ^ 2 t} $, um eine inverse Fourier-Konvertierung durchzuführen. Das Tolle daran ist ** $ x $ und $. Hier wird t $ nicht differenziert. ** Bei der Methode, bei der die Differenzierung wie zuvor verwendet wird, ist es grundlegend, die Differenzierung durch die Differenz zu ersetzen. Daher war es eine absolute Bedingung, dass die Breite der Differenz klein war. Zuerst werden $ x $ und $ f (x) $ unterschieden, aber die Reihenfolge von ** $ \ Delta x $ hat keinen Einfluss auf das Berechnungsergebnis. ** $ t $ wird überhaupt nicht unterschieden Mit anderen Worten, der Fehler sammelt sich nicht an, selbst wenn die Zeit von $ t $ vorverlegt ist. Dies ist ziemlich erstaunlich.

Die letzte Frage ist nun, was $ k $ ist. Welchen Wert nimmt $ k_i $ für $ x_i $, differenziert in ein bestimmtes $ N $? Dies kann durch ernsthafte Betrachtung der diskreten Fourier-Transformation von $ x_i $ gesehen werden.

k_i = 
\begin{cases}
2\pi\frac{i}{N}\hspace{0.5cm} (0 \le i < N/2)\\
2\pi\frac{N-i}{N}\hspace{0.5cm} (N/2 \le i < N)
\end{cases}

Dies ist etwas verwirrend. ** Aber SciPy hat sogar eine Funktion, um dieses $ k_i $ zu generieren. ** Es ist sehr vorsichtig.

Implementierung

Wie ich schon lange gesagt habe, ist die Codierung einfach. ** Die Fourier-Konvertierung ist in scipy.fftpack verfügbar. **

Es gibt zwei Punkte zu beachten.

from scipy.fftpack import fft, ifft, fftfreq

def f(x):
    return np.exp(-x**2)

# set system
N, L = 256, 40.0
x = np.linspace(-L/2, L/2, N)

# set initial function
gaussian = f(x)
plt.plot(x, gaussian, label="t=0")

# set k_i
k = 2 * np.pi * fftfreq(N, d=1/N)/L
time = [2, 4, 8]

for t in time:
    expK = np.exp(-k**2 * t)
    result = ifft(fft(gaussian) * expK)

Es gibt eine for-Anweisung, aber sie ist niedlich ... Diese Art von Zeitschritt ist unvermeidlich. Dieses Mal habe ich sie nur dreimal gedreht. Ich möchte jedoch vorsichtig sein und die Zeit nacheinander wie ** 2s → 4s → 8s vorverlegen Dies bedeutet, dass jeder von ihnen die zeitliche Entwicklung unabhängig berechnet. ** Mit anderen Worten, wenn Sie einen Graphen von 8s möchten, können Sie einfach "t = 8" einsetzen. Dies ist die Differenzmethode. Ist der größte Unterschied.

figure6.png

Sie haben das gleiche Diagramm wie beim letzten Mal erhalten. Diese Methode kann angewendet werden, auch wenn sie sehr leistungsfähig ist und Potenzial hat. Sie wird sich mithilfe der Traberzerlegung zur Symplectic-Methode entwickeln ... aber sie ist zu spezialisiert. Machen wir weiter so.

Nichtlineare Differentialgleichung (Evolution)

Dieser Inhalt ist ein wenig technisch. Wenn Sie nicht interessiert sind, können Sie ihn durchgehen, aber er enthält auch eine kleine wichtige Geschichte.

In der Physik treten häufig nichtlineare Differentialgleichungen auf. Eine solche nichtlineare Gleichung [^ 2]

\left(-\frac{1}{2}\frac{d^2}{dx^2} + g|\psi(x)|\right)\psi(x) = \mu\psi(x)

Dies scheint schwierig zu lösen, wie es früher war, weil ** es $ \ psi (x) $ in der Matrix enthält, wenn es differenziert wird **.

\left[
\frac{1}{2\Delta x^2}
	\begin{pmatrix}
	-2&1&0&\cdots&0\\ 
	1&-2 &1&&0\\
	0 & 1&-2&& \vdots\\
	\vdots&&&\ddots&1\\
	0& \cdots& 0 &1& -2
\end{pmatrix}
+g
\begin{pmatrix}
	|\psi_0|^2&0&0&\cdots&0\\ 
	0&|\psi_1|^2 &0&&0\\
	0 & 0&|\psi_2|^2&& \vdots\\
	\vdots&&&\ddots&0\\
	0& \cdots& 0 &0& |\psi_n|^2
\end{pmatrix}
\right]
\begin{pmatrix}
\psi_0\\
\psi_1\\
\vdots\\
\psi_{n-1}\\
\psi_n
\end{pmatrix}
= T(\psi)\psi
=\mu\psi

Dies bedeutet, dass es nicht mit einem linearen Agonisten wie $ T \ psi $ geschrieben werden kann, was eine Interpretation des Wortes "nichtlinear" ist. Obwohl es oben zwangsweise als Matrix geschrieben wurde, ist $ T (\ psi) Es wird \ psi $. Es ist eine Struktur, die Sie $ \ psi $ suchen möchten, aber die Matrix, die erforderlich ist, um es zu finden, enthält $ \ psi $. Dies wird als selbstkonsistent bezeichnet. tat.

Stellen wir uns nun die obige Differentialgleichung als Eigenwertgleichung mit $ \ mu $ als Eigenwert vor:

T(\psi)\psi = \mu_n \psi

Da es sich um eine selbstkonsistente Gleichung handelt, fallen folgende Ideen ein:

Und es gibt eine Soliton-Lösung in der obigen nichtlinearen Gleichung, und unter ihnen gibt es, wenn $ g <0 $ ist, eine Gauß-ähnliche [^ 3] spezielle Lösung namens Bright Soliton. Also $ \ psi $ Lassen Sie uns diese Idee umsetzen, indem wir Gauß als Ausgangsfunktion von wählen.

Implementierung

def h(x):
    return x * np.exp(-x**2)

# set system
L, N, g = 30, 200, -5
x, dx = np.linspace(-L/2, L/2, N), L / N

# set K matrix
K = np.eye(N, N)
K_sub = np.vstack((K[1:], np.array([0] * N)))
K = (dx**-2 * (2 * K - K_sub - K_sub.T)) * 0.5

# set initial parameter
v, w, mu = np.array([h(x)]).reshape(N, 1), np.array([3]), 0

# self-consistent loop
while abs(w[np.argmin(w)] - mu) > 1e-5:

    # update mu
    mu = w[np.argmin(w)]
    
    # set G matrix
    G = g * np.diag(np.abs(v.T[np.argmin(w)])**2)
    
    # solve
    T = K + G
    w, v = np.linalg.eigh(T)

Diesmal habe ich $ vT [0] $ als zu wiederholende Funktion gewählt. In dieser Differentialgleichung scheint 1 Soliton der Eigenfunktion mit der niedrigsten Energie zu entsprechen [^ 4]. Ich habe ungefähr 16 Mal mit dieser Art von System geloopt. Lassen Sie uns Folgendes zeichnen:

figure7.png

Es sieht aus wie Soliton, aber ich mache mir Sorgen, wenn es nicht Gauß ist. Versuchen Sie es mit Solitons allgemeiner Lösung.

f_s(x) = \frac{b}{a\cosh(x)}

Passen wir zu Gauß:

from scipy.optimize import curve_fit

def bright_soliton(x, a, b):
    return (b / np.cosh(a * x))**2

def gauss_function(x, a, b):
    return a * np.exp(-b * x**2)
	
param_soliton = curve_fit(bright_soliton, x, v.T[0]**2)[0]
param_gauss = curve_fit(gauss_function, x, result)[0]

figure8.png

Das "Ergebnis" und "Bright_Soliton" stimmten genau überein. Es passte nicht zu Gauß, also ist es wahrscheinlich ein Soliton.

Eine etwas schwierige Geschichte wurde fortgesetzt, aber diesmal ist es nicht Soliton usw., sondern in Code wie einer selbstkonsistenten Gleichung, die die Matrix initialisiert und die Berechnung wiederholt, während der Wert geändert wird, Schleife zur Initialisierung nacheinander Dies bedeutet, dass die Verwendung viel Zeit in Anspruch nimmt. ** Es ist schwierig, eine selbstkonsistente Schleife + Matrixinitialisierung durchzuführen. 2 Schleife = 3 Schleife. NumPy ist dafür unverzichtbar. Im vorherigen Basic "Initialisieren Sie jedes Mal eine große Matrix, wenn der Parameter geändert und berechnet wird. Es ist sehr leistungsfähig bei Aufgaben, mit denen Sie fortfahren können ", sagte eine Berechnung wie diese.

Endlich / persönlich

Ich habe über die Implementierung grundlegender numerischer Berechnungen gesprochen. Ich denke, ich werde etwas spezialisierter sein, also höre ich hier auf. Danach hoffe ich, dass Sie die Referenz von NumPy / SciPy gemäß der Berechnung sehen können, die Sie durchführen möchten. Ich denke, ich habe jedoch ein wenig darüber gelassen, so dass ich es in einem anderen Artikel zusammenfassen kann.

Die Motivation für das Schreiben dieses Artikels war in erster Linie, dass ich das Gefühl hatte, dass es nicht viele Informationen über diese Art der Beschleunigung im Internet gibt. Ich folgte dem folgenden Weg.

  1. Beeindruckt von den einfachen Spezifikationen von Python und der Fülle an numerischen Berechnungsbibliotheken. Ich möchte es für meine eigenen Recherchen verwenden, aber es ist sehr langsam, wenn ich den C ++ - Code so ablege, wie er ist. Kann ich mit harten numerischen Berechnungen umgehen? Und fang an zu suchen.

  2. Die Geschichte "NumPy ist schnell!" Ist überall, aber Leute, die von C migrieren, verwenden am Ende die for / while-Anweisung. ** Die Torheit, ein NumPy-Array an die for-Anweisung zu übergeben Ich dachte: "Es ist überhaupt nicht schnell!" ** Es ist lächerlich, wenn ich jetzt darüber nachdenke.

  3. "Denn Anweisung ist langsam!" Wird auch in verschiedenen Bereichen erwähnt, aber ** Personen, die von C migriert sind, wissen nicht, wie man numerischen Berechnungscode schreibt, ohne for / while-Anweisung zu verwenden. ** Es ist wahr, dass die Matrixoperation schnell ist, aber ich war immer noch besorgt darüber, wie man die Matrix vorbereitet, wie man andere als die Matrix berechnet und so weiter.

  4. Es wird gesagt, dass "das Einfügen von Listen für Sätze schneller als gewöhnlich ist!", Aber ehrlich gesagt war es ein kleiner Unterschied.

  5. Und das Ziel war ** Cython, Boost / Python, Numba usw. [^ 5]. ** Diese Tools sind sicherlich schneller, haben aber mehr Einschränkungen, ** immer mehr C-ähnlichen Code Ich werde fortsetzen. **

  6. Das Ergebnis ist "Ist C ++ nicht gut?". Wenn Sie jedoch an Python gewöhnt sind, werden Sie von den esoterischen Sprachspezifikationen von C ++ angewidert sein.

  7. Kehren Sie zu 5 zurück.

Ich wiederhole das seit fast einem Jahr. Ich kann verstehen, wie jede Funktion verwendet wird, indem ich das NumPy-Handbuch lese, aber ich kann nicht viele Erklärungen für numerische Algorithmen finden, die die Verwendung für Anweisungen ausdrücklich gründlich vermeiden. Ich denke das. ** "Verwenden Sie nicht so viel für Aussagen wie möglich! Verwenden Sie die universelle Funktion von NumPy! Ich zeige Ihnen ein Beispiel!" ** Wenn solche Informationen vor Ihnen rollen, können Sie einen solchen Umweg machen. Ich glaube nicht, dass es da war.

Vor diesem Hintergrund beabsichtige ich, es mit dem Ziel zu schreiben, einen kochbuchartigen grundlegenden numerischen Berechnungsalgorithmus zu erstellen. Ich hoffe, er hilft Menschen, die sich Sorgen über die Ausführungsgeschwindigkeit numerischer Berechnungen machen.

[^ 1]: Um genau zu sein, ist die Bedingung, dass "$ f $ analytisch in Fourier konvertiert werden kann" anstelle von "$ f $ kann analytisch integriert werden".

[^ 2]: Es heißt die Gross-Pitaevskii-Gleichung.

[^ 3]: Ähnlich, aber mit völlig anderen Eigenschaften als Gauß.

[^ 4]: Es besteht eine gute Chance, dass es sich je nach zu lösendem Modell ändert, und möglicherweise sind Versuche und Irrtümer erforderlich.

[^ 5]: Ich hoffe, ich kann diese auch eines Tages kommentieren.

Recommended Posts

Beschleunigung der numerischen Berechnung mit NumPy / SciPy: Anwendung 2
Beschleunigung der numerischen Berechnung mit NumPy / SciPy: Anwendung 1
Beschleunigung der numerischen Berechnung mit NumPy / SciPy: Aufnehmen gefallener Ohren
Beschleunigung numerischer Berechnungen mit NumPy: Basics
Sehen Sie, wie schnell Sie mit NumPy / SciPy beschleunigen können
[Python] Beschleunigung der Verarbeitung mit Cache-Tools
[Wissenschaftlich-technische Berechnung mit Python] Lösen (verallgemeinerter) Eigenwertprobleme mit numpy / scipy mithilfe von Bibliotheken
Verwenden von Intel MKL mit NumPy / SciPy (Version November 2019)
[Competition Pro] Beschleunigung der DP mithilfe der kumulierten Summe
Versuchen Sie es mit scipy