En mathématiques , et plus particulièrement en analyse , l'exponentielle d'une matrice est une fonction généralisant la fonction exponentielle aux matrices et aux endomorphismes par le calcul fonctionnel . Elle fait en particulier le pont entre un groupe de Lie et son algèbre de Lie .
Théorème et définition — La série de matrices de terme général
M
n
(
C
)
⟶
M
n
(
C
)
A
⟼
1
k
!
A
k
{\displaystyle {\begin{array}{ccc}\mathrm {M} _{n}(\mathbb {C} )&\longrightarrow &\mathrm {M} _{n}(\mathbb {C} )\\A&\longmapsto &{\frac {1}{k!}}A^{k}\end{array}}}
converge normalement sur toute partie bornée de
M
n
(
C
)
{\displaystyle \mathrm {M} _{n}(\mathbb {C} )}
[ 1] .
On appelle alors exponentielle l'application
exp
:
M
n
(
C
)
→
M
n
(
C
)
,
A
↦
e
A
=
∑
k
∈
N
1
k
!
A
k
{\displaystyle \exp :\mathrm {M} _{n}(\mathbb {C} )\to \mathrm {M} _{n}(\mathbb {C} ),\quad A\mapsto \mathrm {e} ^{A}=\sum _{k\in \mathbb {N} }{\frac {1}{k!}}A^{k}}
.
Pour n = 1, on retrouve la définition de l'exponentielle complexe .
Sauf indication contraire, X , Y , etc. désignent des matrices n × n complexes (à coefficients complexes).
L'exponentielle de la matrice nulle est la matrice identité :
e
0
=
I
{\displaystyle \mathrm {e} ^{0}=I}
;
Le déterminant de l'exponentielle d'une matrice est égal à l'exponentielle de sa trace :
det
(
e
X
)
=
e
tr
(
X
)
{\displaystyle \det \left(\mathrm {e} ^{X}\right)=\mathrm {e} ^{\operatorname {tr} (X)}}
;
si Y est une matrice inversible , alors
e
Y
X
Y
−
1
=
Y
e
X
Y
−
1
{\displaystyle \mathrm {e} ^{YXY^{-1}}=Y\,\mathrm {e} ^{X}Y^{-1}}
;
l'exponentielle de matrice vérifie la limite :
e
X
=
lim
k
→
∞
(
I
+
X
k
)
k
{\displaystyle \mathrm {e} ^{X}=\lim _{k\to \infty }\left(I+{\frac {X}{k}}\right)^{\!k}}
;
e
(
X
+
Y
)
=
lim
n
→
∞
(
e
X
n
e
Y
n
)
n
{\displaystyle \mathrm {e} ^{(X+Y)}=\lim _{n\to \infty }\left(\mathrm {e} ^{\frac {X}{n}}\mathrm {e} ^{\frac {Y}{n}}\right)^{\!n}}
(formule de Trotter-Kato ) ;
il existe un polynôme d'endomorphisme PX (dépendant de X ) tel que
e
X
=
P
X
(
X
)
{\displaystyle \mathrm {e} ^{X}=P_{X}(X)}
.
L'exponentielle induit une surjection continue de
A
n
(
R
)
{\displaystyle \mathrm {A_{n}} (\mathbb {R} )}
(ensemble des matrices antisymétriques ) dans
S
O
n
(
R
)
{\displaystyle \mathrm {SO_{n}} (\mathbb {R} )}
(groupe spécial orthogonal des matrices orthogonales ) de déterminant
1
{\displaystyle 1}
[ 2] .
Elle induit de plus un homéomorphisme de
S
n
(
R
)
{\displaystyle \mathrm {S_{n}} (\mathbb {R} )}
(ensemble des matrices symétriques) dans
S
n
+
+
(
R
)
{\displaystyle \mathrm {{S_{n}}^{++}} (\mathbb {R} )}
(ensemble des matrices symétriques définies positives ).
La transposée , la conjuguée et l'adjointe d'une matrice X sont notées
X
T
{\displaystyle X^{\mathsf {T}}}
,
X
¯
{\displaystyle {\overline {X}}}
et
X
†
{\displaystyle X^{\dagger }}
.
L'exponentielle de la transposée d'une matrice est la transposée de l'exponentielle de la matrice :
e
X
T
=
(
e
X
)
T
{\displaystyle \mathrm {e} ^{X^{\mathsf {T}}}=\left(\mathrm {e} ^{X}\right)^{\!{\mathsf {T}}}}
. Il s'ensuit que :
si X est symétrique (
X
T
=
X
{\displaystyle X^{\mathsf {T}}=X}
), alors eX l'est aussi :
(
e
X
)
T
=
e
X
{\displaystyle \left(\mathrm {e} ^{X}\right)^{\!{\mathsf {T}}}=\mathrm {e} ^{X}}
;
si X est antisymétrique (
X
T
=
−
X
{\displaystyle X^{\mathsf {T}}=-X}
) et réelle (à coefficients réels), alors eX est orthogonale :
(
e
X
)
−
1
=
(
e
X
)
T
{\displaystyle \left(\mathrm {e} ^{X}\right)^{-1}=\left(\mathrm {e} ^{X}\right)^{\mathsf {T}}}
.
L'exponentielle de la conjuguée d'une matrice est la conjuguée de l'exponentielle de la matrice :
e
X
¯
=
e
X
¯
{\displaystyle \mathrm {e} ^{\overline {X}}={\overline {\mathrm {e} ^{X}}}}
et donc, compte tenu de la propriété précédente :
L'exponentielle de l'adjointe d'une matrice est l'adjointe de l'exponentielle de la matrice :
e
X
†
=
(
e
X
)
†
{\displaystyle \mathrm {e} ^{X^{\dagger }}=\left(\mathrm {e} ^{X}\right)^{\!\dagger }}
. Il s'ensuit que :
si X est hermitienne (
X
†
=
X
{\displaystyle X^{\dagger }=X}
), alors eX l'est aussi :
(
e
X
)
†
=
e
X
{\displaystyle \left(\mathrm {e} ^{X}\right)^{\!\dagger }=\mathrm {e} ^{X}}
;
si X est antihermitienne (
X
†
=
−
X
{\displaystyle X^{\dagger }=-X}
), alors eX est unitaire :
(
e
X
)
−
1
=
(
e
X
)
†
{\displaystyle \left(\mathrm {e} ^{X}\right)^{-1}=\left(\mathrm {e} ^{X}\right)^{\dagger }}
.
Le commutateur de X et Y est noté [X , Y ] (= XY – YX ) .
Si [X , Y ] = 0 (les matrices commutent) alors
e
X
e
Y
=
e
X
+
Y
(
=
e
Y
e
X
)
{\displaystyle \mathrm {e} ^{X}\mathrm {e} ^{Y}=\mathrm {e} ^{X+Y}\;(=\mathrm {e} ^{Y}\mathrm {e} ^{X})}
.
Plus généralement, en supposant seulement que [X , Y ] commute avec X et Y ,
e
X
e
Y
=
e
X
+
Y
+
1
2
[
X
,
Y
]
{\displaystyle \mathrm {e} ^{X}\mathrm {e} ^{Y}=\mathrm {e} ^{X+Y+{\frac {1}{2}}[X,Y]}}
(formule de Glauber ).
Soit
Φ
(
t
)
=
e
X
t
e
Y
t
{\displaystyle \Phi (t)=\mathrm {e} ^{Xt}\mathrm {e} ^{Yt}}
. Alors,
Φ
(
0
)
=
I
{\displaystyle \Phi (0)=I}
et
Φ
′
(
t
)
=
X
e
X
t
e
Y
t
+
e
X
t
Y
e
Y
t
=
(
X
+
e
X
t
Y
e
−
X
t
)
Φ
(
t
)
=
(
X
+
Y
+
[
e
X
t
,
Y
]
e
−
X
t
)
Φ
(
t
)
{\displaystyle \Phi '(t)=X\mathrm {e} ^{Xt}\mathrm {e} ^{Yt}+\mathrm {e} ^{Xt}Y\mathrm {e} ^{Yt}=(X+\mathrm {e} ^{Xt}Y\mathrm {e} ^{-Xt})\Phi (t)=\left(X+Y+\left[\mathrm {e} ^{Xt},Y\right]\mathrm {e} ^{-Xt}\right)\Phi (t)}
.
Remarquons d'abord la relation
[
A
,
F
(
B
)
]
=
[
A
,
B
]
F
′
(
B
)
{\displaystyle [A,F(B)]=[A,B]\,F'(B)}
, pour toutes matrices A et B telles que B commute avec [A , B ] et pour toute fonction F développable en série entière[ 5] . En effet, pour tout
n
∈
N
{\displaystyle n\in \mathbb {N} }
, la formule générale
[
A
,
B
n
]
=
∑
0
≤
k
<
n
B
k
[
A
,
B
]
B
n
−
1
−
k
{\displaystyle [A,B^{n}]=\sum _{0\leq k<n}B^{k}\left[A,B\right]B^{n-1-k}}
devient ici :
[
A
,
B
n
]
=
n
[
A
,
B
]
B
n
−
1
{\displaystyle [A,B^{n}]=n[A,B]B^{n-1}}
.
Appliquons cette relation au commutateur
[
Y
,
e
X
t
]
{\displaystyle \left[Y,\mathrm {e} ^{Xt}\right]}
:
[
Y
,
e
X
t
]
=
t
[
Y
,
X
]
e
X
t
{\displaystyle \left[Y,\mathrm {e} ^{Xt}\right]=t[Y,X]\mathrm {e} ^{Xt}}
, ou encore
[
e
X
t
,
Y
]
=
t
[
X
,
Y
]
e
X
t
{\displaystyle \left[\mathrm {e} ^{Xt},Y\right]=t[X,Y]\mathrm {e} ^{Xt}}
, donc
Φ
′
(
t
)
=
(
X
+
Y
+
t
[
X
,
Y
]
)
Φ
(
t
)
{\displaystyle \Phi '(t)=\left(X+Y+t[X,Y]\right)\Phi (t)}
.
On intègre en utilisant que [X , Y ] commute avec X + Y :
Φ
(
t
)
=
Φ
(
0
)
e
t
(
X
+
Y
)
+
t
2
2
[
X
,
Y
]
=
e
t
(
X
+
Y
)
+
t
2
2
[
X
,
Y
]
{\displaystyle \Phi (t)=\Phi (0)\mathrm {e} ^{t(X+Y)+{\frac {t^{2}}{2}}[X,Y]}=\mathrm {e} ^{t(X+Y)+{\frac {t^{2}}{2}}[X,Y]}}
En prenant t = 1 :
e
X
e
Y
=
e
X
+
Y
+
1
2
[
X
,
Y
]
{\displaystyle \mathrm {e} ^{X}\mathrm {e} ^{Y}=\mathrm {e} ^{X+Y+{\frac {1}{2}}[X,Y]}}
.
Encore plus généralement, la formule de Baker-Campbell-Hausdorff donne l'expression de
e
X
+
Y
−
e
X
e
Y
{\displaystyle {\rm {e}}^{X+Y}-\mathrm {e} ^{X}\mathrm {e} ^{Y}}
, plus précisément d'un logarithme de eX eY , par une série ne faisant intervenir que X , Y et leurs commutateurs. Les premiers termes sont[ 6] , [ 7] :
log
(
e
X
e
Y
)
=
X
+
Y
+
1
2
[
X
,
Y
]
+
1
12
[
X
,
[
X
,
Y
]
]
+
1
12
[
Y
,
[
Y
,
X
]
]
−
1
24
[
X
,
[
Y
,
[
X
,
Y
]
]
]
+
…
{\displaystyle \log(\mathrm {e} ^{X}\mathrm {e} ^{Y})=X+Y+{\frac {1}{2}}[X,Y]+{\frac {1}{12}}[X,[X,Y]]+{\frac {1}{12}}[Y,[Y,X]]-{\frac {1}{24}}[X,[Y,[X,Y]]]+\dots }
La réciproque est cependant fausse : pour
X
=
(
0
0
0
2
i
π
)
,
Y
=
(
0
0
1
2
i
π
)
,
{\displaystyle X={\begin{pmatrix}0&0\\0&2\mathrm {i} \pi \end{pmatrix}},Y={\begin{pmatrix}0&0\\1&2\mathrm {i} \pi \end{pmatrix}},}
on a
e
X
=
e
Y
=
e
X
+
Y
=
(
1
0
0
1
)
{\displaystyle \mathrm {e} ^{X}=\mathrm {e} ^{Y}=\mathrm {e} ^{X+Y}={\begin{pmatrix}1&0\\0&1\end{pmatrix}}}
mais ces deux matrices ne commutent pas.
L'exponentielle d'une matrice est toujours inversible . L'inverse de eX est donné par e-X . Cette fonction est donc une application
exp
:
M
n
(
C
)
→
G
L
n
(
C
)
{\displaystyle \exp :\mathrm {M} _{n}(\mathbb {C} )\to \mathrm {GL} _{n}(\mathbb {C} )}
de l'ensemble des matrices n × n vers le groupe général linéaire , c'est-à-dire le groupe de toutes les matrices inversibles. Cette application est surjective .
Pour deux matrices X et Y , nous avons :
‖
e
X
+
Y
−
e
X
‖
≤
‖
Y
‖
e
‖
X
‖
e
‖
Y
‖
,
{\displaystyle \|\mathrm {e} ^{X+Y}\!-\mathrm {e} ^{X}\|\leq \|Y\|\,\mathrm {e} ^{\|X\|}\,\mathrm {e} ^{\|Y\|},}
où || · || désigne une norme matricielle arbitraire. Il suit que l'application exponentielle est continue et lipschitzienne sur tout sous-ensemble compact de
M
n
(
C
)
{\displaystyle \mathrm {M} _{n}(\mathbb {C} )}
.
L'application
exp
:
M
n
(
C
)
→
G
L
n
(
C
)
{\displaystyle \exp :\mathrm {M} _{n}(\mathbb {C} )\to \mathrm {GL} _{n}(\mathbb {C} )}
est même de classe
C
∞
{\displaystyle C^{\infty }}
.
Sa différentielle en 0 est l'identité et elle réalise un difféomorphisme entre un voisinage de 0 et un voisinage de l'identité.
L'application :
t
↦
e
t
X
,
t
∈
R
{\displaystyle t\mapsto {\rm {e}}^{tX},\quad t\in \mathbb {R} }
définit une courbe de classe
C
∞
{\displaystyle C^{\infty }}
dans le groupe linéaire qui passe par l'identité en t = 0. Cette courbe est en fait un sous-groupe de Lie commutatif à un paramètre de
G
L
n
(
C
)
{\displaystyle \mathrm {GL} _{n}(\mathbb {C} )}
puisque :
e
t
1
X
e
t
2
X
=
e
(
t
1
+
t
2
)
X
=
e
t
2
X
e
t
1
X
{\displaystyle \mathrm {e} ^{t_{1}X}\mathrm {e} ^{t_{2}X}=\mathrm {e} ^{(t_{1}+t_{2})X}=\mathrm {e} ^{t_{2}X}\mathrm {e} ^{t_{1}X}}
.
La dérivée de cette courbe au point t est donnée par :
d
d
t
e
t
X
=
X
e
t
X
=
e
t
X
X
(
1
)
{\displaystyle {\frac {\mathrm {d} }{\mathrm {d} t}}\mathrm {e} ^{tX}=X\,\mathrm {e} ^{tX}=\mathrm {e} ^{tX}X\qquad (1)}
(la dérivée au point t = 0 est la matrice X , ce qui revient à dire que X engendre ce sous-groupe à un paramètre)
En effet, plus généralement, la différentielle de l'application exponentielle en une matrice X est donnée par :
exp
X
′
(
H
)
=
∑
i
,
j
∈
N
X
i
H
X
j
(
i
+
j
+
1
)
!
=
∑
i
,
j
∈
N
B
(
i
+
1
,
j
+
1
)
X
i
i
!
H
X
j
j
!
=
∫
0
1
e
(
1
−
α
)
X
H
e
α
X
d
α
{\displaystyle \exp '_{X}(H)=\sum _{i,j\in \mathbb {N} }{\frac {X^{i}HX^{j}}{(i+j+1)!}}=\sum _{i,j\in \mathbb {N} }\mathrm {B} (i+1,j+1){\frac {X^{i}}{i!}}H{\frac {X^{j}}{j!}}=\int _{0}^{1}{\rm {e}}^{(1-\alpha )X}H{\rm {e}}^{\alpha X}{\rm {d}}\alpha }
où B désigne la fonction bêta , d'où :
si
X
H
=
H
X
,
exp
X
′
(
H
)
=
H
e
X
{\displaystyle {\text{si }}XH=HX,\quad \exp '_{X}(H)=H\,\mathrm {e} ^{X}}
.
Dans le plan euclidien muni d'un repère orthonormé , considérons la matrice de rotation d'angle π / 2 :
R
⊥
=
(
0
−
1
1
0
)
{\displaystyle R_{\perp }={\begin{pmatrix}0&-1\\1&0\end{pmatrix}}}
.
Alors :
e
θ
R
⊥
=
(
cos
θ
−
sin
θ
sin
θ
cos
θ
)
=
R
(
θ
)
{\displaystyle \mathrm {e} ^{\theta \,R_{\perp }}={\begin{pmatrix}\cos \theta &-\sin \theta \\\sin \theta &\cos \theta \end{pmatrix}}=R(\theta )}
est la matrice de rotation d'angle θ [ 8] .
Soit la matrice :
P
G
=
(
0
1
0
0
)
{\displaystyle P_{\mathrm {G} }={\begin{pmatrix}0&1\\0&0\end{pmatrix}}}
.
Alors :
e
−
v
P
G
=
(
1
−
v
0
1
)
=
T
G
(
v
)
{\displaystyle \mathrm {e} ^{-v\,P_{\mathrm {G} }}={\begin{pmatrix}1&-v\\0&1\end{pmatrix}}=T_{\mathrm {G} }(v)}
est la matrice de transformation de Galilée dans le plan (x , t ) pour un déplacement de vitesse v sur l'axe Ox [ 8] : x' = x – vt , t' = t .
Soit la matrice :
P
L
=
(
0
1
1
0
)
{\displaystyle P_{\mathrm {L} }={\begin{pmatrix}0&1\\1&0\end{pmatrix}}}
.
Alors :
e
−
φ
P
L
=
(
cosh
φ
−
sinh
φ
−
sinh
φ
cosh
φ
)
=
T
L
(
φ
)
{\displaystyle \mathrm {e} ^{-\varphi P_{\mathrm {L} }}={\begin{pmatrix}\cosh \varphi &-\sinh \varphi \\-\sinh \varphi &\cosh \varphi \end{pmatrix}}=T_{\mathrm {L} }(\varphi )}
est la matrice de transformation de Lorentz dans le plan (x , ct ) pour un déplacement de rapidité φ sur l'axe Ox [ 8] .
Soit
n
^
{\displaystyle {\hat {n}}}
un vecteur unitaire de cosinus directeurs α , β , γ (
n
^
=
α
e
→
x
+
β
e
→
y
+
γ
e
→
z
{\displaystyle {\hat {n}}=\alpha \,{\vec {e}}_{x}+\beta \,{\vec {e}}_{y}+\gamma \,{\vec {e}}_{z}\,}
avec α 2 + β 2 + γ 2 = 1 ), et soit la matrice[ a] :
N
n
^
=
(
0
−
γ
β
γ
0
−
α
−
β
α
0
)
{\displaystyle N_{\hat {n}}={\begin{pmatrix}0&-\gamma &\beta \\\gamma &0&-\alpha \\-\beta &\alpha &0\end{pmatrix}}}
.
Alors :
e
θ
N
n
^
=
R
n
^
(
θ
)
{\displaystyle \mathrm {e} ^{\theta \,N_{\hat {n}}}=R_{\hat {n}}(\theta )}
est la matrice de rotation d'angle θ autour d'un axe Δ de vecteur unitaire
n
^
{\displaystyle {\hat {n}}}
[ 8] .
En géologie structurale , on s'intéresse à la déformation finie résultant, au bout d'un certain temps, d'une déformation progressive[ 9] :
r
→
f
=
D
⋅
r
→
0
{\displaystyle {\vec {r}}_{\mathrm {f} }=D\cdot {\vec {r}}_{0}}
,
d
r
→
d
t
=
L
(
t
)
⋅
r
→
(
t
)
{\displaystyle {\frac {\mathrm {d} {\vec {r}}}{\mathrm {d} t}}=L(t)\cdot {\vec {r}}(t)}
,
où
r
→
(
t
)
=
O
M
→
{\displaystyle {\vec {r}}(t)={\overrightarrow {\mathrm {OM} }}}
désigne le vecteur position par rapport à un point matériel arbitraire choisi comme origine (qui peut suivre n'importe quelle trajectoire entre les instants t 0 et t f ),
r
→
0
{\displaystyle {\vec {r}}_{0}}
la position initiale (à
t
=
t
0
{\displaystyle t=t_{0}}
) et
r
→
f
{\displaystyle {\vec {r}}_{\mathrm {f} }}
la position finale (à t = t f ). D est la « matrice de déformation finie » et L (t ) la « matrice de déformation progressive ».
Si L (t ) est une matrice constante L , alors :
D
=
e
(
t
−
t
0
)
L
{\displaystyle D=\mathrm {e} ^{(t-t_{0})L}}
.
Si L (t ) varie mais commute avec sa dérivée
(
d
L
/
d
t
)
{\displaystyle ({\mathrm {d} L}/{\mathrm {d} t})}
[ b] , alors la formule précédente se généralise en[ 9] , [ c] :
D
=
e
∫
t
0
t
f
L
(
τ
)
d
τ
{\displaystyle D=\mathrm {e} ^{\int _{t_{0}}^{t_{\mathrm {f} }}L(\tau )\,\mathrm {d} \tau }}
.
Dans le cas général on ne sait exprimer D que sous la forme d'un développement en série [ 10] , [ d] .
Le calcul d'une exponentielle de matrice n'est pas a priori un problème facile. Cependant, dans certains cas, et notamment ceux d'une matrice diagonale et d'une matrice nilpotente, il ne présente aucune difficulté. Une fois cette remarque faite, le cas général peut se traiter en se ramenant aux deux cas précédents.
Si D est une matrice diagonale , c'est-à-dire :
D
=
(
d
1
0
…
0
0
d
2
…
0
⋮
⋮
⋱
⋮
0
0
…
d
n
)
{\displaystyle D={\begin{pmatrix}d_{1}&0&\ldots &0\\0&d_{2}&\ldots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\ldots &d_{n}\end{pmatrix}}}
,
alors son exponentielle est obtenue en calculant l'exponentielle de chacun des termes de la diagonale principale :
e
D
=
(
e
d
1
0
…
0
0
e
d
2
…
0
⋮
⋮
⋱
⋮
0
0
…
e
d
n
)
{\displaystyle \mathrm {e} ^{D}={\begin{pmatrix}\mathrm {e} ^{d_{1}}&0&\ldots &0\\0&\mathrm {e} ^{d_{2}}&\ldots &0\\\vdots &\vdots &\ddots &\vdots \\0&0&\ldots &\mathrm {e} ^{d_{n}}\end{pmatrix}}}
.
Si A est une matrice diagonalisable , c'est-à-dire :
A
=
P
D
P
−
1
{\displaystyle A=PDP^{-1}}
où D est diagonale, alors
e
A
=
P
e
D
P
−
1
{\displaystyle \mathrm {e} ^{A}=P\mathrm {e} ^{D}P^{-1}}
.
L'application exponentielle préserve ainsi les espaces propres , soit les sous-espaces engendrés par les vecteurs colonnes de P .
De plus, les valeurs propres de eA sont les exponentielles de celles de A , soit les éléments diagonaux de eD .
Une matrice N est nilpotente si N q = 0 pour un entier q . Dans ce cas, son exponentielle eN se calcule directement à partir de son développement en série, puisque celui-ci ne comporte alors qu'un nombre fini de termes :
e
N
=
I
+
N
+
N
2
2
!
+
N
3
3
!
+
⋯
+
N
q
−
1
(
q
−
1
)
!
{\displaystyle \operatorname {e} ^{N}=I+N+{\frac {N^{2}}{2!}}+{\frac {N^{3}}{3!}}+\cdots +{\frac {N^{q-1}}{(q-1)!}}}
.
Lorsque le polynôme minimal d'une matrice X est scindé (ce qui est en particulier toujours le cas pour les matrices à coefficients complexes), la décomposition de Dunford donne
X
=
A
+
N
{\displaystyle X=A+N\,}
où
A est diagonalisable ;
N est nilpotente ;
A commute avec N .
Dès lors, le calcul de l'exponentielle de X se réduit aux deux cas précédents :
e
X
=
e
A
+
N
=
e
A
e
N
{\displaystyle \mathrm {e} ^{X}=\mathrm {e} ^{A+N}=\mathrm {e} ^{A}\mathrm {e} ^{N}}
.
On peut aussi faire appel à la réduction de Jordan : soit J la forme de Jordan de X , et P la matrice de passage . Alors,
e
X
=
P
e
J
P
−
1
{\displaystyle \mathrm {e} ^{X}=P\mathrm {e} ^{J}P^{-1}}
.
Puisque
J
=
J
a
1
(
λ
1
)
⊕
J
a
2
(
λ
2
)
⊕
⋯
⊕
J
a
n
(
λ
n
)
{\displaystyle J=J_{a_{1}}(\lambda _{1})\oplus J_{a_{2}}(\lambda _{2})\oplus \cdots \oplus J_{a_{n}}(\lambda _{n})}
,
e
J
{\displaystyle \mathrm {e} ^{J}}
=
exp
(
J
a
1
(
λ
1
)
⊕
J
a
2
(
λ
2
)
⊕
⋯
⊕
J
a
n
(
λ
n
)
)
{\displaystyle =\exp \left(J_{a_{1}}(\lambda _{1})\oplus J_{a_{2}}(\lambda _{2})\oplus \cdots \oplus J_{a_{n}}(\lambda _{n})\right)}
=
exp
(
J
a
1
(
λ
1
)
)
⊕
exp
(
J
a
2
(
λ
2
)
)
⊕
⋯
⊕
exp
(
J
a
k
(
λ
k
)
)
{\displaystyle =\exp \left(J_{a_{1}}(\lambda _{1})\right)\oplus \exp \left(J_{a_{2}}(\lambda _{2})\right)\oplus \cdots \oplus \exp \left(J_{a_{k}}(\lambda _{k})\right)}
.
En conséquence, il faut seulement connaître la méthode pour calculer l'exponentielle d'un bloc de Jordan. Chacun est de la forme
J
a
(
λ
)
=
λ
I
+
N
{\displaystyle J_{a}(\lambda )=\lambda I+N}
où N est une matrice nilpotente. L'exponentielle du bloc est donnée par
e
λ
I
+
N
=
e
λ
e
N
{\displaystyle \mathrm {e} ^{\lambda I+N}=\mathrm {e} ^{\lambda }\mathrm {e} ^{N}}
.
Soit la matrice
B
=
(
21
17
6
−
5
−
1
−
6
4
4
16
)
{\displaystyle B={\begin{pmatrix}21&17&6\\-5&-1&-6\\4&4&16\end{pmatrix}}}
qui a la forme de Jordan
J
=
P
−
1
B
P
=
(
4
0
0
0
16
1
0
0
16
)
{\displaystyle J=P^{-1}BP={\begin{pmatrix}4&0&0\\0&16&1\\0&0&16\end{pmatrix}}}
et la matrice de passage
P
=
(
−
1
4
2
5
4
1
4
−
2
−
1
4
0
4
0
)
{\displaystyle P={\begin{pmatrix}-{\frac {1}{4}}&2&{\frac {5}{4}}\\{\frac {1}{4}}&-2&-{\frac {1}{4}}\\0&4&0\end{pmatrix}}}
,
d'inverse
P
−
1
=
(
1
5
2
0
0
1
4
1
1
0
)
{\displaystyle P^{-1}={\begin{pmatrix}1&5&2\\0&0&{\frac {1}{4}}\\1&1&0\end{pmatrix}}}
.
Maintenant,
J
=
J
1
(
4
)
⊕
J
2
(
16
)
et
e
B
=
P
e
J
P
−
1
=
P
(
e
J
1
(
4
)
⊕
e
J
2
(
16
)
)
P
−
1
{\displaystyle J=J_{1}(4)\oplus J_{2}(16)\qquad {\text{et}}\qquad \mathrm {e} ^{B}=P\mathrm {e} ^{J}P^{-1}=P\left(\mathrm {e} ^{J_{1}(4)}\oplus \mathrm {e} ^{J_{2}(16)}\right)P^{-1}}
.
L'exponentielle de la matrice 1×1 J 1 (4) = [4] est simplement la matrice 1×1 [e4 ].
L'exponentielle de la matrice 2×2 J 2 (16) peut se calculer par la formule eλI +N = eλ eN mentionnée ci-dessus ; on obtient
exp
(
16
I
+
(
0
1
0
0
)
)
=
e
16
(
(
1
0
0
1
)
+
(
0
1
0
0
)
+
1
2
!
(
0
0
0
0
)
+
⋯
)
=
(
e
16
e
16
0
e
16
)
{\displaystyle \exp \left(16I+{\begin{pmatrix}0&1\\0&0\end{pmatrix}}\right)=\mathrm {e} ^{16}\left({\begin{pmatrix}1&0\\0&1\end{pmatrix}}+{\begin{pmatrix}0&1\\0&0\end{pmatrix}}+{1 \over 2!}{\begin{pmatrix}0&0\\0&0\end{pmatrix}}+\cdots \right)={\begin{pmatrix}\mathrm {e} ^{16}&\mathrm {e} ^{16}\\0&\mathrm {e} ^{16}\end{pmatrix}}}
,
d'où
e
B
=
P
(
e
4
0
0
0
e
16
e
16
0
0
e
16
)
P
−
1
=
1
4
(
13
e
16
−
e
4
13
e
16
−
5
e
4
2
e
16
−
2
e
4
−
9
e
16
+
e
4
−
9
e
16
+
5
e
4
−
2
e
16
+
2
e
4
16
e
16
16
e
16
4
e
16
)
{\displaystyle \mathrm {e} ^{B}=P{\begin{pmatrix}\mathrm {e} ^{4}&0&0\\0&\mathrm {e} ^{16}&\mathrm {e} ^{16}\\0&0&\mathrm {e} ^{16}\end{pmatrix}}P^{-1}={1 \over 4}{\begin{pmatrix}13\mathrm {e} ^{16}-\mathrm {e} ^{4}&13\mathrm {e} ^{16}-5\mathrm {e} ^{4}&2\mathrm {e} ^{16}-2\mathrm {e} ^{4}\\-9\mathrm {e} ^{16}+\mathrm {e} ^{4}&-9\mathrm {e} ^{16}+5\mathrm {e} ^{4}&-2\mathrm {e} ^{16}+2\mathrm {e} ^{4}\\16\mathrm {e} ^{16}&16\mathrm {e} ^{16}&4\mathrm {e} ^{16}\end{pmatrix}}}
.
Une des premières applications de l'exponentielle de matrices est la résolution des équations différentielles ordinaires . En effet, de l'équation ci-dessus, on déduit que la solution de :
d
d
t
y
(
t
)
=
A
y
(
t
)
,
y
(
0
)
=
y
0
{\displaystyle {\frac {\mathrm {d} }{\mathrm {d} t}}y(t)=Ay(t),\quad y(0)=y_{0}}
,
où A est une matrice, est donnée par
y
(
t
)
=
e
A
t
y
0
{\displaystyle y(t)=\mathrm {e} ^{At}y_{0}}
.
L'exponentielle d'une matrice peut aussi servir à résoudre les équations non homogènes :
d
d
t
y
(
t
)
=
A
y
(
t
)
+
b
(
t
)
,
y
(
0
)
=
y
0
{\displaystyle {\frac {\mathrm {d} }{\mathrm {d} t}}y(t)=Ay(t)+b(t),\quad y(0)=y_{0}}
.
En multipliant par e−At , nous avons
e
−
A
t
d
d
t
y
(
t
)
−
e
−
A
t
A
y
(
t
)
=
e
−
A
t
b
{\displaystyle \mathrm {e} ^{-At}{\frac {\mathrm {d} }{\mathrm {d} t}}y(t)-\mathrm {e} ^{-At}Ay(t)=\mathrm {e} ^{-At}b}
d
d
t
(
e
−
A
t
y
)
=
e
−
A
t
b
{\displaystyle {\frac {\mathrm {d} }{\mathrm {d} t}}(\mathrm {e} ^{-At}y)=\mathrm {e} ^{-At}b}
.
La résolution du système se ramène donc au calcul de eAt .
Il n'existe pas de solution explicite pour les équations différentielles de la forme :
d
d
t
y
(
t
)
=
A
(
t
)
y
(
t
)
,
y
(
0
)
=
y
0
{\displaystyle {\frac {\mathrm {d} }{\mathrm {d} t}}y(t)=A(t)\,y(t),\quad y(0)=y_{0}}
où A n'est pas constant, mais le développement de Magnus (en) donne la solution sous la forme d'une somme infinie.
Soit le système
{
x
′
=
2
x
−
y
+
z
y
′
=
3
y
−
z
z
′
=
2
x
+
y
+
3
z
.
{\displaystyle \left\{{\begin{matrix}x'&=&2x&-y&+z\\y'&=&&3y&-z\\z'&=&2x&+y&+3z.\end{matrix}}\right.}
La matrice associée est
M
=
(
2
−
1
1
0
3
−
1
2
1
3
)
{\displaystyle M={\begin{pmatrix}2&-1&1\\0&3&-1\\2&1&3\end{pmatrix}}}
et son exponentielle est
e
t
M
=
(
−
1
1
1
1
0
−
1
1
−
1
1
)
(
e
t
J
2
(
2
)
⊕
e
t
J
1
(
4
)
)
(
1
/
2
1
1
/
2
1
1
0
1
/
2
0
1
/
2
)
=
(
(
e
2
t
(
1
−
2
t
)
+
e
4
t
)
/
2
−
t
e
2
t
(
−
e
2
t
+
e
4
t
)
/
2
(
e
2
t
(
1
+
2
t
)
−
e
4
t
)
/
2
e
2
t
(
1
+
t
)
(
e
2
t
−
e
4
t
)
/
2
(
e
2
t
(
−
1
+
2
t
)
+
e
4
t
)
/
2
t
e
2
t
(
e
2
t
+
e
4
t
)
/
2
)
.
{\displaystyle {\begin{aligned}\mathrm {e} ^{tM}&={\begin{pmatrix}-1&1&1\\1&0&-1\\1&-1&1\end{pmatrix}}\left(\mathrm {e} ^{tJ_{2}(2)}\oplus \mathrm {e} ^{tJ_{1}(4)}\right){\begin{pmatrix}1/2&1&1/2\\1&1&0\\1/2&0&1/2\end{pmatrix}}\\&={\begin{pmatrix}(\mathrm {e} ^{2t}(1-2t)+\mathrm {e} ^{4t})/2&-t\mathrm {e} ^{2t}&(-\mathrm {e} ^{2t}+\mathrm {e} ^{4t})/2\\(\mathrm {e} ^{2t}(1+2t)-\mathrm {e} ^{4t})/2&\mathrm {e} ^{2t}(1+t)&(\mathrm {e} ^{2t}-\mathrm {e} ^{4t})/2\\(\mathrm {e} ^{2t}(-1+2t)+\mathrm {e} ^{4t})/2&t\mathrm {e} ^{2t}&(\mathrm {e} ^{2t}+\mathrm {e} ^{4t})/2\end{pmatrix}}.\end{aligned}}}
La solution générale du système est donc
(
x
y
z
)
=
C
1
2
(
e
2
t
(
1
−
2
t
)
+
e
4
t
e
2
t
(
1
+
2
t
)
−
e
4
t
e
2
t
(
−
1
+
2
t
)
+
e
4
t
)
+
C
2
(
−
t
e
2
t
e
2
t
(
1
+
t
)
t
e
2
t
)
+
C
3
2
(
−
e
2
t
+
e
4
t
e
2
t
−
e
4
t
e
2
t
+
e
4
t
)
{\displaystyle {\begin{pmatrix}x\\y\\z\end{pmatrix}}={\frac {C_{1}}{2}}{\begin{pmatrix}\mathrm {e} ^{2t}(1-2t)+\mathrm {e} ^{4t}\\\mathrm {e} ^{2t}(1+2t)-\mathrm {e} ^{4t}\\\mathrm {e} ^{2t}(-1+2t)+\mathrm {e} ^{4t}\end{pmatrix}}+C_{2}{\begin{pmatrix}-t\mathrm {e} ^{2t}\\\mathrm {e} ^{2t}(1+t)\\t\mathrm {e} ^{2t}\end{pmatrix}}+{\frac {C_{3}}{2}}{\begin{pmatrix}-\mathrm {e} ^{2t}+\mathrm {e} ^{4t}\\\mathrm {e} ^{2t}-\mathrm {e} ^{4t}\\\mathrm {e} ^{2t}+\mathrm {e} ^{4t}\end{pmatrix}}}
c'est-à-dire, en posant
a
=
C
1
/
2
{\displaystyle a=C_{1}/2}
,
b
=
C
2
{\displaystyle b=C_{2}}
et
c
=
C
3
/
2
{\displaystyle c=C_{3}/2}
:
x
=
a
(
e
2
t
(
1
−
2
t
)
+
e
4
t
)
−
b
t
e
2
t
+
c
(
−
e
2
t
+
e
4
t
)
y
=
a
(
e
2
t
(
1
+
2
t
)
−
e
4
t
)
+
b
e
2
t
(
1
+
t
)
+
c
(
e
2
t
−
e
4
t
)
z
=
a
(
e
2
t
(
−
1
+
2
t
)
+
e
4
t
)
+
b
t
e
2
t
+
c
(
e
2
t
+
e
4
t
)
.
{\displaystyle {\begin{aligned}x&=a(\mathrm {e} ^{2t}(1-2t)+\mathrm {e} ^{4t})-bt\mathrm {e} ^{2t}+c(-\mathrm {e} ^{2t}+\mathrm {e} ^{4t})\\y&=a(\mathrm {e} ^{2t}(1+2t)-\mathrm {e} ^{4t})+b\mathrm {e} ^{2t}(1+t)+c(\mathrm {e} ^{2t}-\mathrm {e} ^{4t})\\z&=a(\mathrm {e} ^{2t}(-1+2t)+\mathrm {e} ^{4t})+bt\mathrm {e} ^{2t}+c(\mathrm {e} ^{2t}+\mathrm {e} ^{4t}).\end{aligned}}}
Pour une équation non homogène, on peut utiliser une méthode semblable à la variation de la constante .
Nous cherchons une solution de la forme y p (t ) = exp(tA )z (t ) :
y
p
′
=
(
e
t
A
)
′
z
(
t
)
+
e
t
A
z
′
(
t
)
=
A
e
t
A
z
(
t
)
+
e
t
A
z
′
(
t
)
=
A
y
p
(
t
)
+
e
t
A
z
′
(
t
)
{\displaystyle \mathbf {y} _{p}'=(\mathrm {e} ^{tA})'\mathbf {z} (t)+\mathrm {e} ^{tA}\mathbf {z} '(t)=A\mathrm {e} ^{tA}\mathbf {z} (t)+\mathrm {e} ^{tA}\mathbf {z} '(t)=A\mathbf {y} _{p}(t)+\mathrm {e} ^{tA}\mathbf {z} '(t)}
Avec y p comme solution :
e
t
A
z
′
(
t
)
=
b
(
t
)
{\displaystyle \mathrm {e} ^{tA}\mathbf {z} '(t)=\mathbf {b} (t)}
z
′
(
t
)
=
(
e
t
A
)
−
1
b
(
t
)
{\displaystyle \mathbf {z} '(t)=(\mathrm {e} ^{tA})^{-1}\mathbf {b} (t)}
z
(
t
)
=
∫
0
t
e
−
u
A
b
(
u
)
d
u
+
c
{\displaystyle \mathbf {z} (t)=\int _{0}^{t}\mathrm {e} ^{-uA}\mathbf {b} (u)\,\mathrm {d} u+\mathbf {c} }
.
Alors,
y
p
=
e
t
A
∫
0
t
e
−
u
A
b
(
u
)
d
u
+
e
t
A
c
=
∫
0
t
e
(
t
−
u
)
A
b
(
u
)
d
u
+
e
t
A
c
{\displaystyle \mathbf {y} _{p}=\mathrm {e} ^{tA}\int _{0}^{t}\mathrm {e} ^{-uA}\mathbf {b} (u)\,\mathrm {d} u+\mathrm {e} ^{tA}\mathbf {c} =\int _{0}^{t}\mathrm {e} ^{(t-u)A}\mathbf {b} (u)\,\mathrm {d} u+\mathrm {e} ^{tA}\mathbf {c} }
où c dépend des conditions initiales.
Soit le système
{
x
′
=
2
x
−
y
+
z
+
e
2
t
y
′
=
3
y
−
1
z
z
′
=
2
x
+
y
+
3
z
+
e
2
t
.
{\displaystyle \left\{{\begin{matrix}x'&=&2x&-y&+z&+\mathrm {e} ^{2t}\\y'&=&&3y&-1z&\\z'&=&2x&+y&+3z&+\mathrm {e} ^{2t}.\end{matrix}}\right.}
Nous avons donc
M
=
(
2
−
1
1
0
3
−
1
2
1
3
)
et
b
=
e
2
t
(
1
0
1
)
{\displaystyle M={\begin{pmatrix}2&-1&1\\0&3&-1\\2&1&3\end{pmatrix}}\qquad {\text{et}}\qquad \mathbf {b} =\mathrm {e} ^{2t}{\begin{pmatrix}1\\0\\1\end{pmatrix}}}
.
Comme auparavant, la somme de la solution homogène et de la solution particulière donne la solution générale. La solution homogène étant connue, il suffit de trouver la solution particulière.
y
p
=
e
t
∫
0
t
e
−
u
A
(
e
2
u
0
e
2
u
)
d
u
+
e
t
A
c
=
e
t
∫
0
t
(
2
e
u
−
2
u
e
2
u
−
2
u
e
2
u
0
−
2
e
u
+
2
(
u
+
1
)
e
2
u
2
(
u
+
1
)
e
2
u
0
2
u
e
2
u
2
u
e
2
u
2
e
u
)
(
e
2
u
0
e
2
u
)
d
u
+
e
t
A
c
=
e
t
∫
0
t
(
e
2
u
(
2
e
u
−
2
u
e
2
u
)
e
2
u
(
−
2
e
u
+
2
(
1
+
u
)
e
2
u
)
2
e
3
u
+
2
u
e
4
u
)
d
u
+
e
t
A
c
=
e
t
(
−
1
24
e
3
t
(
3
e
t
(
4
t
−
1
)
−
16
)
1
24
e
3
t
(
3
e
t
(
4
t
+
4
)
−
16
)
1
24
e
3
t
(
3
e
t
(
4
t
−
1
)
−
16
)
)
+
(
2
e
t
−
2
t
e
2
t
−
2
t
e
2
t
0
−
2
e
t
+
2
(
t
+
1
)
e
2
t
2
(
t
+
1
)
e
2
t
0
2
t
e
2
t
2
t
e
2
t
2
e
t
)
(
c
1
c
2
c
3
)
,
{\displaystyle {\begin{aligned}\mathbf {y} _{p}&=\mathrm {e} ^{t}\int _{0}^{t}\mathrm {e} ^{-uA}{\begin{pmatrix}\mathrm {e} ^{2u}\\0\\\mathrm {e} ^{2u}\end{pmatrix}}\,\mathrm {d} u+e^{tA}\mathbf {c} \\&=\mathrm {e} ^{t}\int _{0}^{t}{\begin{pmatrix}2\mathrm {e} ^{u}-2u\mathrm {e} ^{2u}&-2u\mathrm {e} ^{2u}&0\\\\-2\mathrm {e} ^{u}+2(u+1)\mathrm {e} ^{2u}&2(u+1)\mathrm {e} ^{2u}&0\\\\2u\mathrm {e} ^{2u}&2u\mathrm {e} ^{2u}&2\mathrm {e} ^{u}\end{pmatrix}}{\begin{pmatrix}\mathrm {e} ^{2u}\\0\\\mathrm {e} ^{2u}\end{pmatrix}}\,\mathrm {d} u+\mathrm {e} ^{tA}\mathbf {c} \\&=\mathrm {e} ^{t}\int _{0}^{t}{\begin{pmatrix}\mathrm {e} ^{2u}(2\mathrm {e} ^{u}-2u\mathrm {e} ^{2u})\\\\\mathrm {e} ^{2u}(-2\mathrm {e} ^{u}+2(1+u)\mathrm {e} ^{2u})\\\\2\mathrm {e} ^{3u}+2u\mathrm {e} ^{4u}\end{pmatrix}}\,\mathrm {d} u+\mathrm {e} ^{tA}\mathbf {c} \\&=\mathrm {e} ^{t}{\begin{pmatrix}-{1 \over 24}\mathrm {e} ^{3t}(3\mathrm {e} ^{t}(4t-1)-16)\\\\{1 \over 24}\mathrm {e} ^{3t}(3\mathrm {e} ^{t}(4t+4)-16)\\\\{1 \over 24}\mathrm {e} ^{3t}(3\mathrm {e} ^{t}(4t-1)-16)\end{pmatrix}}+{\begin{pmatrix}2\mathrm {e} ^{t}-2t\mathrm {e} ^{2t}&-2t\mathrm {e} ^{2t}&0\\\\-2\mathrm {e} ^{t}+2(t+1)\mathrm {e} ^{2t}&2(t+1)\mathrm {e} ^{2t}&0\\\\2t\mathrm {e} ^{2t}&2t\mathrm {e} ^{2t}&2\mathrm {e} ^{t}\end{pmatrix}}{\begin{pmatrix}c_{1}\\c_{2}\\c_{3}\end{pmatrix}},\end{aligned}}}
expression qui peut être simplifiée pour obtenir la solution particulière cherchée.
↑ Cette matrice est celle qui permet d'écrire le produit vectoriel comme une application linéaire :
N
n
^
r
→
=
n
^
∧
r
→
.
{\displaystyle N_{\hat {n}}\,{\vec {r}}={\hat {n}}\wedge {\vec {r}}.}
↑ C'est notamment le cas quand L (t ) est proportionnelle à une matrice constante (L (t ) = L 0 f (t ) ), ou bien encore si elle est diagonale .
↑ Pour vérifier que cette expression est bien une (la) solution du système différentiel et des conditions initiales ci-dessus, il suffit de calculer
d
r
→
/
d
t
{\displaystyle {\mathrm {d} {\vec {r}}}/{\mathrm {d} t}}
en appliquant la définition de l'exponentielle d'une matrice :
e
A
=
I
+
A
+
A
2
2
+
⋯
+
A
n
n
!
+
⋯
{\displaystyle \mathrm {e} ^{A}=I+A+{\frac {A^{2}}{2}}+\cdots +{\frac {A^{n}}{n!}}+\cdots }
.
↑ On connaît une solution analytique fermée dans quelques rares cas où
L
(
t
)
{\displaystyle L(t)}
ne commute pas avec sa dérivée, notamment celui d'une matrice triangulaire [ 10] .
↑ Voir par exemple le chapitre « Exponentielle d'une matrice » sur Wikiversité .
↑ Maxime Bouchereau, « Exponentielle de matrices antisymétriques » [PDF]
↑ C. Cohen-Tannoudji , B. Diu et F. Laloë , Mécanique quantique [détail de l’édition ] , vol. 1, p. 174-175 (en anglais sur Google Livres et p. 158-159 en allemand sur Google Livres ).
↑ Albert Messiah , Mécanique quantique [détail des éditions ] , vol. 1, p. 442 de la traduction en anglais .
↑ Cohen-Tannoudji, Diu et Laloë , p. 171-172.
↑ Roger Godement , Introduction à la théorie des groupes de Lie , Berlin, Springer , 1982 , 305 p. (ISBN 3-540-20034-7 , lire en ligne ) , p. 263 .
↑ Pour plus de termes, voir par exemple (en) Janusz Czyż , Paradoxes of Measures and Dimensions Originating in Felix Hausdorff's Ideas , World Scientific, 1994 , 738 p. (ISBN 978-981-02-0189-0 , lire en ligne ) , p. 421 .
↑ a b c et d Jean-Pierre Provost et Gérard Vallée, Les maths en physique : La physique à travers le filtre des mathématiques , Paris, Éditions Dunod , coll. « Sciences Sup », mars 2004 , 1re éd. , 331 p. (ISBN 2-10-004652-7 ) , p. 101-102 .
↑ a et b (en) Ariel Provost, Cécile Buisson et Olivier Merle, « From progressive to finite deformation and back », Journal of Geophysical Research: Solid Earth , vol. 109, no B2, février 2004 , p. 1-11 , article no B02405 (DOI 10.1029/2001JB001734 , lire en ligne , consulté le 10 juin 2018 ) .
↑ a et b Daniel Pham, Techniques du Calcul Matriciel , Paris, Dunod, 1962 , 387 p. , p. 232-235 .
Sur les autres projets Wikimedia :
(en) Roger A. Horn et Charles R. Johnson , Topics in Matrix Analysis , Cambridge University Press , 1991 , 607 p. (ISBN 0-521-46713-6 , lire en ligne )
Xavier Merlin, Algèbre , Ellipses , coll. « Methodix », 2007 , 400 p. (ISBN 978-2-7298-9555-6 et 2-7298-9555-8 )
(en) Cleve Moler et Charles Van Loan , « Nineteen Dubious Ways to Compute the Exponential of a Matrix », SIAM Review , vol. 20, no 4, 1978 (DOI 10.1137/1020098 )
(en) Cleve Moler et Charles Van Loan , « Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later », SIAM Review , vol. 45, no 1, 2003 (DOI 10.1137/S00361445024180 )
(en) Roger B. Sidje, « Expokit: a software package for computing matrix exponentials », ACM TOMS , vol. 24, no 1, mars 1998 (DOI 10.1145/285861.285868 ) — Code source