introduction

Pour les opérateurs non commutables $ X et Y $, la formule de Lie-Trotter suivante contient [^ nom].

\exp(h(X+Y)) = \left( \exp(hX/n)\exp(hY/n) \right)^n + O\left(\frac{h^2}{n}\right)

Ici, $ h $ est le nombre c et est souvent utilisé comme pas de temps dans le calcul numérique, donc ce qui suit est appelé le pas de temps. $ n $ est le nombre de décomposition. Le but de cet article est de confirmer correctement cette erreur de coupure.

La source est https://github.com/kaityo256/lie-trotter-sample À.

Ajouté le 14 juillet 2017: La source a été légèrement modifiée afin que vous puissiez également essayer la décomposition symétrique quadratique. Pour plus de détails, reportez-vous à l'article suivant Décomposition symétrique du second ordre dans la formule de Lie-Trotter.

Explication de ce que fait le script

Considérons une matrice carrée $ X $, $ Y $ de la dimension $ d $ appropriée. Faisons une matrice symétrique pour que les valeurs propres soient réelles. Calculez strictement $ \ exp (h (X + Y)) $, $ \ exp (hX) $, $ \ exp (hY) $, etc., et calculez l'erreur en utilisant la formule de Lie-Trotter. Pour l'erreur, j'utiliserai la norme Frobenius.

Faire une vraie matrice symétrique

Tout d'abord, créez au hasard une matrice symétrique de dimension $ d $ appropriée.

import numpy as np
import math
from scipy import linalg

d = 2
np.random.seed(1)
x1 = np.random.rand(d,d)
x2 = np.random.rand(d,d)
X = x1.dot(x1.T)
Y = x2.dot(x2.T)
Z = X + Y

Créez une matrice aléatoire appropriée et définissez le produit de la translocation sur $ X $ ou $ Y $. De cette manière, la valeur propre devient réelle car elle devient une matrice symétrique. Supposons que la somme de $ X $ et $ Y $ soit $ Z $.

Calcul de la fonction exponentielle de la matrice

Créez une fonction qui place la matrice sur l'épaule de l'exposant à l'intervalle de temps spécifié $ h $. Pour ce faire, commencez par diagonaliser.

X = U_x \Lambda_x U_x^{-1}

Parce qu'il est facile de mettre une matrice diagonale sur l'épaule d'une fonction exponentielle

\exp(h X) = U_x \exp(h \Lambda_x) U_x^{-1}

Peut être calculé comme. L'écriture en Python ressemble à ceci.

def calc_exp(X, h):
    rx, Ux = linalg.eigh(X)
    Ux_inv = linalg.inv(Ux)
    tx = np.diag(np.array([math.exp(h * i) for i in rx]))
    eX = Ux.dot(tx.dot(Ux_inv))
    return eX

Expliquez ce que vous faites

Appelez linalg.eigh car il se nourrit d'une matrice symétrique
Mettre le vecteur de valeur propre sur l'épaule de la fonction exponentielle dans la notation d'inclusion de la liste
Faites-en une matrice avec numpy.diag
Annuler avec $ U_x $ ou $ U_x ^ {-1} $

Se sentir comme ça.

Nombre spécifié de décompositions Décomposition au trotteur

Ce sera facile si la fonction exponentielle de la matrice peut être calculée dans l'intervalle de temps spécifié. Par exemple, dans le cas de $ h $ en incréments de temps et $ n $ en nombre de décomposition, créez d'abord $ \ exp (h X / n) $ et $ \ exp (h Y / n) $. Tout ce que vous avez à faire est de multiplier la matrice unitaire $ n $ fois.

    eX = calc_exp(X,h/n)
    eY = calc_exp(Y,h/n)
    S = np.diag(np.ones(d))
    eXeY = eX.dot(eY)
    for i in range(n):
        S = S.dot(eXeY)

Comparez la matrice approximative $ S $ ainsi créée avec l'évaluation exacte $ \ exp (h Z) $. Comparons ici avec Frobenius Norm. En réunissant tout ce qui précède, la fonction ressemble à ceci.

def trotter(X,Y,Z,h,n):
    eZ = calc_exp(Z,h)
    eX = calc_exp(X,h/n)
    eY = calc_exp(Y,h/n)
    S = np.diag(np.ones(d))
    eXeY = eX.dot(eY)
    for i in range(n):
        S = S.dot(eXeY)
    return linalg.norm(eZ - S)/linalg.norm(eZ)

résultat

Dans le cas bidimensionnel, changeons $ h $ et $ n $ pour voir ce qui arrive à l'erreur de coupure.

Tout d'abord, la dépendance $ h $ de l'erreur de troncature. Soit $ n $ 1,2,4, et évaluez l'erreur de coupure avec différentes valeurs $ h $ pour chacune.

On peut voir que la dépendance $ h $ de l'erreur est $ (h ^ 2) $ quelle que soit la valeur de $ n $. Un $ n $ plus grand réduit l'erreur, qui est proportionnelle à 1 $ / n $. Fixez-le à $ h = 1.0 $ et évaluez l'erreur de coupure pour divers $ n $.

On voit que l'erreur est certainement $ O (1 / n) $.

Résumé

J'ai évalué l'erreur de coupure de la décomposition du trotteur. Jusqu'à présent, le cas $ \ exp (h (A + B)) \ sim \ exp (hA) \ exp (hB) $ de $ n = 1 $ était appelé "décomposition du premier ordre", donc l'erreur est en quelque sorte $ O. Je pensais que c'était (h) $, mais maintenant je sais que c'était $ O (h ^ 2) $ quelle que soit la valeur de $ n $.

[^ name]: Certains oncles effrayants peuvent venir dire la bonne chose à propos de ce nom officiel. Pour le moment, si l'opérateur est borné, la formule de Lie est utilisée, si elle n'est pas bornée, la formule Trotter est utilisée, et l'algorithme d'intégration de chemin utilisant cela s'appelle la décomposition de Suzuki-Trotter. J'appelle souvent cela simplement le démontage du trotteur, mais ...

[PYTHON] Examiner l'erreur de coupure de la formule de Lie-Trotter