[PYTHON] Accélération du calcul numérique avec NumPy / SciPy: Application 2

Cible

Il est écrit pour les débutants de Python et NumPy. Il peut être particulièrement utile pour ceux qui visent le calcul numérique de systèmes physiques correspondant à "Je peux utiliser le langage C mais j'ai récemment commencé Python" ou "Je ne comprends pas comment écrire comme Python". Peut être.

De plus, il peut y avoir des descriptions incorrectes en raison de mon manque d'étude. Veuillez me pardonner.

Résumé

Accélération du calcul numérique avec NumPy: Basics Accélération du calcul numérique à l'aide de NumPy / SciPy: Application 1

Je vais suivre la méthode de calcul numérique de base, mais cette fois, elle comprend également un peu de contenu avancé.

Algèbre / Équation transcendantale

Les équations algébriques sont ce qu'on appelle des équations solubles à la main. Les équations transcendantales sont des noms assez exagérés, mais ce sont des équations qui ne peuvent pas être résolues par des méthodes algébriques.

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

Ceci est un enfant. Cette équation dit: "Vous avez besoin de $ \ frac {x} {2} $ pour connaître $ \ sin (x) $, et de $ pour connaître $ \ frac {x} {2} $. On l'appelle aussi une «équation auto-cohérente» parce qu'elle a une structure que «\ sin (x) $ est nécessaire». Bien sûr, elle peut être résolue numériquement même si elle ne peut pas être résolue algébriquement. Les algorithmes peuvent être "méthode Newton" ou "dichotomie". Je vais omettre les détails de l'algorithme, mais ce n'est pas si difficile, donc si vous ne le connaissez pas, veuillez le vérifier. Ces implémentations sont simples, mais au final, substitution séquentielle Il est inévitable d'utiliser une boucle for car c'est comme faire. La méthode de Newton converge essentiellement rapidement, mais une fonction avec de bonnes propriétés est nécessaire à condition que ce soit une fonction différenciable. D'autre part, la dichotomie peut toujours être trouvée s'il y a une solution dans l'intervalle spécifié, mais le nombre d'itérations jusqu'à la convergence est grande.

Cependant, comme vous l'avez peut-être déjà remarqué, ** ceux-ci sont déjà implémentés dans un module appelé scipy.optimize. ** Les boucles For sont complètes dans la fonction et s'exécutent rapidement. ..

la mise en oeuvre

Découvrez où se trouve approximativement la solution de l'équation ci-dessus:

figure5.png

Vous pouvez voir qu'il s'agit de $ x \ in (1.5, 2.0) $. Je résoudrai ceci dans une dichotomie:

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

C'est très simple. S'il y a deux ou plusieurs solutions dans l'intervalle, vous ne savez pas laquelle va converger, alors assurez-vous d'en avoir une. Bien sûr, vous pouvez également utiliser des équations algébriques.

Conversion de Fourier

En parlant de conversion de Fourier, nous sommes familiers avec le traitement du signal, etc., mais je vais continuer cette fois aussi l'histoire de la physique. Veuillez comprendre qu'il y aura un peu plus de discussion mathématique. Rappelez-vous l'équation de diffusion traitée dans l'article précédent. S'il te plait donne moi:

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

Cet enfant pourrait être résolu avec la transformée de Fourier. Voyons comment. Nous définissons la transformée de Fourier et la transformée de Fourier inverse comme suit:

\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)]

Remplaçons cette définition par l'équation de diffusion, qui nous permet d'effectuer une différenciation $ x $:

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

De plus, si vous supprimez l'intégration $ k $, cela devient juste une équation différentielle du premier ordre pour $ t $, donc elle peut être facilement résolue:

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

Et je vais faire une conversion de Fourier inverse de ceci:

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

Je ne sais pas s'il reste à $ \ tilde {f} $, donc je vais le corriger à $ f $:

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

Si l'équation de la condition initiale $ f (x, 0) $ est connue et qu'elle peut être intégrée analytiquement, alors la solution analytique [^ 1] peut être obtenue, si elle ne peut pas être intégrée analytiquement. Après tout, ce n'est pas résolu numériquement.

Maintenant, une simple description de ce flux est la suivante:

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

Après avoir converti $ f (x, 0) $ en Fourier, multipliez par $ e ^ {-k ^ 2 t} $ pour effectuer une conversion de Fourier inverse. L'avantage à ce sujet est ** $ x $ et $. C'est là que t $ n'est pas différencié. ** Dans la méthode utilisant la différenciation comme précédemment, il est fondamental de remplacer la différenciation par la différence, donc c'était une condition absolue que la largeur de la différence soit petite. Différencie d'abord $ x $ et $ f (x) $, mais l'ordre de ** $ \ Delta x $ n'affecte pas le résultat du calcul. ** $ t $ n'est pas différencié en premier lieu En d'autres termes, l'erreur ne s'accumule pas même si l'heure de $ t $ est avancée, ce qui est assez étonnant.

Maintenant, la dernière question est de savoir ce qu'est $ k $. Quelle valeur prend $ k_i $ pour $ x_i $ différencié en un certain $ N $? Cela peut être vu en considérant sérieusement la transformation de Fourier discrète de $ x_i $,

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}

C'est un peu déroutant. ** Mais SciPy a même une fonction pour générer ce $ k_i $. ** C'est très prudent.

la mise en oeuvre

Comme je l'ai dit depuis longtemps, le codage est simple. ** La conversion de Fourier est disponible dans scipy.fftpack. **

Il y a deux points à noter.

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)

Il y a une instruction for, mais c'est mignon ... Ce genre de pas de temps est inévitable. Cette fois je ne l'ai tourné que 3 fois. Cependant, je veux faire attention, avancer le temps séquentiellement comme ** 2s → 4s → 8s Cela signifie que chacun d'eux calcule l'évolution du temps indépendamment. ** En d'autres termes, si vous voulez un graphique de 8s, vous pouvez simplement remplacer «t = 8». C'est la méthode de la différence. C'est la plus grande différence.

figure6.png

Vous avez le même graphe que la dernière fois. Cette méthode peut être appliquée même lorsqu'elle est très puissante et a du potentiel, et elle évoluera vers la méthode Symplectic utilisant la décomposition du trotteur ... mais elle est trop spécialisée. Continuons jusqu'à ici.

Équation différentielle non linéaire (évolution)

Ce contenu est un peu technique, si vous n'êtes pas intéressé, vous pouvez le parcourir, mais il comprend également une petite histoire importante.

Les équations différentielles non linéaires apparaissent souvent en physique. Une de ces équations non linéaires [^ 2]

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

Cela semble difficile à résoudre de manière traditionnelle, car ** contient $ \ psi (x) $ dans la matrice lorsqu'il est différencié **.

\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

Cela signifie qu'il ne peut pas être écrit avec un agoniste linéaire tel que $ T \ psi $, qui est une interprétation du mot "non linéaire". Bien qu'il soit écrit de force ci-dessus sous forme de matrice, $ T (\ psi) Il devient \ psi $. C'est une structure que vous voulez trouver $ \ psi $, mais la matrice requise pour la trouver contient $ \ psi $. Cela s'appelle auto-cohérent. fait.

Maintenant, considérons l'équation différentielle ci-dessus comme une équation de valeur propre avec $ \ mu $ comme valeur propre:

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

Puisqu'il s'agit d'une équation auto-cohérente, les idées suivantes viennent à l'esprit:

Et il y a une solution de Soliton dans l'équation non linéaire ci-dessus, et parmi elles, quand $ g <0 $, il y a une solution spéciale de type gaussien [^ 3] appelée Soliton brillant. Donc, $ \ psi $ Implémentons cette idée en choisissant Gaussien comme fonction initiale de.

la mise en oeuvre

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)

Cette fois j'ai choisi $ vT [0] $ comme fonction à répéter. Dans cette équation différentielle, il semble que 1 Soliton correspond à la fonction propre d'énergie la plus basse [^ 4]. La condition de convergence est imposée à $ \ mu $ J'ai bouclé environ 16 fois avec ce type de système. Tracons ceci:

figure7.png

Cela ressemble à Soliton, mais je suis inquiet si ce n'est pas gaussien. Essayez la solution générale de Soliton.

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

Faisons correspondre avec Gaussian:

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

Le "result" et "bright_soliton" correspondaient exactement. Il ne correspondait pas à la Gaussienne, donc c'est probablement un Soliton.

Une petite histoire difficile a continué, mais cette fois ce n'est pas Soliton etc., mais dans un code tel qu'une équation auto-cohérente où la matrice est initialisée et le calcul est répété en changeant la valeur **, il boucle à l'initialisation un par un. Cela signifie que cela prendra beaucoup de temps si vous utilisez. ** Il est difficile de faire quelque chose comme boucle auto-cohérente + initialisation de la matrice 2 boucle = 3 boucle. NumPy est indispensable pour le faire. Dans le précédent Basic, "Initialisez une énorme matrice à chaque fois que le paramètre est modifié et calculez. Il est très puissant dans les tâches que vous pouvez exécuter », a-t-il déclaré.

Enfin / personnel

J'ai parlé de l'implémentation des calculs numériques de base. Je pense que je serai un peu plus spécialisé, donc c'est tout. Après ça, j'espère que vous pourrez voir la référence NumPy / SciPy en fonction du calcul que vous voulez faire. Je pense, mais j'en ai laissé un peu de côté, alors je peux le résumer dans un autre article.

La motivation pour écrire cet article en premier lieu était que je sentais qu'il n'y avait pas beaucoup d'informations sur ce genre d'accélération sur le net.J'ai suivi le chemin suivant.

  1. Impressionné par les spécifications simples de Python et l'abondance de bibliothèques de calculs numériques. Je veux l'utiliser pour mes propres recherches, mais il est très lent si je laisse tomber le code C ++ tel quel. Puis-je gérer des calculs numériques difficiles? Et commencez à chercher.

  2. L'histoire "NumPy est rapide!" Est partout, mais les personnes qui migrent de C finissent par utiliser l'instruction for / while. ** La folie de passer un tableau NumPy à l'instruction for , J'ai pensé: "Ce n'est pas du tout rapide!" ** C'est ridicule quand j'y pense maintenant.

  3. "For statement is slow!" On parle également dans divers domaines, mais ** Les personnes qui ont migré de C ne savent pas comment écrire un code de calcul numérique sans utiliser l'instruction for / while. ** Il est vrai que le fonctionnement de la matrice est rapide, mais je m'inquiétais toujours de savoir comment préparer la matrice, comment calculer autre chose que la matrice, etc.

  4. On dit que "l'inclusion de liste est plus rapide que l'ordinaire pour les phrases!", Mais honnêtement, c'était une légère différence.

  5. Et la destination était ** Cython, Boost / Python, Numba, etc. [^ 5]. ** Ces outils sont certes plus rapides mais ont plus de restrictions, ** de plus en plus de code C-like. Je vais continuer. **

  6. Le résultat est "C ++ n'est-il pas bon?" Cependant, si vous êtes habitué à Python, vous serez dégoûté par les spécifications du langage ésotérique du C ++.

  7. Revenez à 5.

Je répète cela depuis presque un an. Je peux comprendre comment utiliser chaque fonction en lisant le manuel NumPy, mais je ne trouve pas beaucoup d'explications sur l'algorithme numérique qui déclare explicitement que l'instruction for n'est pas utilisée. Je pense que. ** "Ne l'utilisez pas autant que possible pour les déclarations! Utilisez la fonction universelle de NumPy! Je vais vous montrer un exemple!" ** Si des informations comme celles-ci roulent devant vous, vous pouvez faire un tel détour. Je ne pense pas que c'était là.

Dans ce contexte, j'ai l'intention de l'écrire dans le but de créer un algorithme de calcul numérique de base semblable à un livre de cuisine et j'espère qu'il aidera les gens qui s'inquiètent de la vitesse d'exécution des calculs numériques.

[^ 1]: Pour être exact, la condition est que "$ f $ peut être converti analytiquement en Fourier" au lieu de "$ f $ peut être intégré analytiquement".

[^ 2]: Cela s'appelle l'équation de Gross-Pitaevskii.

[^ 3]: Similaire, mais avec des propriétés complètement différentes de la Gaussienne.

[^ 4]: Il y a de fortes chances que cela change en fonction du modèle à résoudre, et des essais et erreurs peuvent être nécessaires.

[^ 5]: J'espère que je pourrai aussi commenter cela un jour.

Recommended Posts

Accélération du calcul numérique avec NumPy / SciPy: Application 2
Accélération du calcul numérique avec NumPy / SciPy: Application 1
Accélération du calcul numérique avec NumPy / SciPy: ramasser les oreilles tombées
Accélérer les calculs numériques avec NumPy: principes de base
Découvrez la puissance de l'accélération avec NumPy / SciPy
[Python] Accélération du traitement à l'aide des outils de cache
[Calcul scientifique / technique par Python] Résolution de problèmes de valeurs propres (généralisés) en utilisant numpy / scipy, en utilisant des bibliothèques
Utilisation d'Intel MKL avec NumPy / SciPy (version de novembre 2019)
[Competition Pro] Accélération du DP en utilisant la somme cumulée
Essayez d'utiliser scipy