§E20.18行列関数の数値計算

最終更新

定数係数の線形微分方程式系y′=Ayy'=Ay、y(0)=y0y(0)=y_0の解は、行列指数関数を用いてy(t)=etAy0y(t)=e^{tA}y_0と表される。実正方行列AAが正則行列VVと実対角行列Λ\LambdaによってA=VΛV−1A=V\Lambda V^{-1}と表されるならばeA=VeΛV−1e^A=Ve^\Lambda V^{-1}であり、eΛe^\LambdaはΛ\Lambdaの対角成分の指数関数を並べた対角行列である。しかし、数値計算で得られる固有値と固有ベクトルは近似であり、近似V^\widehat V、Λ^\widehat\Lambdaから作ったV^eΛ^V^−1\widehat Ve^{\widehat\Lambda}\widehat V^{-1}はeAe^Aに一致するとは限らない。固有ベクトルを並べた行列の条件が悪いと、固有値の小さな誤差が、計算値では、AAに同程度の大きさの摂動を加えたときのeAe^Aの変化よりはるかに大きな誤差として現れることがある。

対角化を経ずにeAe^Aを近似する方法に、s∈N≥0s\in\Nに対してAAを2−sA2^{-s}Aに縮小し、そこで指数関数の Taylor 多項式TmT_mを評価して、得た行列をss回二乗するものがある。これを引数の縮小と繰り返す二乗という。eAe^Aはe2−sAe^{2^{-s}A}の2s2^s乗に等しく、縮小した引数の作用素ノルムは2−s∥A∥2^{-s}\lVert A\rVertであるから、A≠0A\ne0ならば、Taylor 近似の打切り誤差の上界はssを大きくすると小さくなる。他方、打切り誤差と二乗に入る丸めの誤差は、二乗を繰り返すたびに伝わる。

この方法の誤差は、厳密算術では∥A∥\lVert A\rVert、次数mm、縮小回数ssで表される上界をもつ。Taylor 多項式の代わりに[1/1][1/1]型の Padé 近似を用いる場合にも、縮小した引数の作用素ノルムが22未満であれば同じ形の上界が得られ、これらの上界は、y(t)y(t)の計算値の相対誤差の上界へ移される。本記事では、行列指数関数の数値計算の方法とその誤差の評価について解説する。

1 行列指数関数と対角化による計算

補題 1.1.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、定義域と値域にこのノルムを入れたX∈Mn(R)X\in M_n(\R)の作用素ノルムを∥X∥\lVert X\rVertと書く。任意のX,Y∈Mn(R)X,Y\in M_n(\R)に対して、∥I∥=1\lVert I\rVert=1、∥XY∥≤∥X∥∥Y∥\lVert XY\rVert\le\lVert X\rVert\lVert Y\rVert、∥eX∥≤e∥X∥\lVert e^X\rVert\le e^{\lVert X\rVert}が成り立つ。

証明.h∈Rnh\in\R^nをとる。Ih=hIh=hであるから∥I∥=1\lVert I\rVert=1である。§E20.2 補題 1.2 (1)により∥XYh∥≤∥X∥∥Yh∥≤∥X∥∥Y∥∥h∥\lVert XYh\rVert\le\lVert X\rVert\lVert Yh\rVert\le\lVert X\rVert\lVert Y\rVert\lVert h\rVertであり、∥h∥=1\lVert h\rVert=1を満たすhhにわたる上限をとって∥XY∥≤∥X∥∥Y∥\lVert XY\rVert\le\lVert X\rVert\lVert Y\rVertを得る。したがって作用素ノルムは劣乗法的であり、§E10.10 命題 1.2により∥eX∥≤∥I∥+e∥X∥−1=e∥X∥\lVert e^X\rVert\le\lVert I\rVert+e^{\lVert X\rVert}-1=e^{\lVert X\rVert}である。▨

補題 1.2.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。任意のX,E∈Mn(R)X,E\in M_n(\R)に対して

∥eX+E−eX∥≤∥E∥e∥X∥+∥E∥\lVert e^{X+E}-e^X\rVert\le\lVert E\rVert e^{\lVert X\rVert+\lVert E\rVert}

が成り立つ。

証明. 演習とする（問題 7.1）。▨

命題 1.3.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。A∈Mn(R)A\in M_n(\R)とし、V^∈Mn(R)\widehat V\in M_n(\R)を正則行列、Λ^=diag⁡(μ1,…,μn)\widehat\Lambda=\operatorname{diag}(\mu_1,\dots,\mu_n)を実対角行列とする。R:=AV^−V^Λ^R:=A\widehat V-\widehat V\widehat\Lambda、E:=−RV^−1E:=-R\widehat V^{-1}と置く。

  1. V^eΛ^V^−1=eA+E\widehat Ve^{\widehat\Lambda}\widehat V^{-1}=e^{A+E}であり、eΛ^=diag⁡(eμ1,…,eμn)e^{\widehat\Lambda}=\operatorname{diag}(e^{\mu_1},\dots,e^{\mu_n})である。特にA=V^Λ^V^−1A=\widehat V\widehat\Lambda\widehat V^{-1}ならばeA=V^eΛ^V^−1e^A=\widehat Ve^{\widehat\Lambda}\widehat V^{-1}である。
  2. ∥E∥≤∥R∥∥V^−1∥\lVert E\rVert\le\lVert R\rVert\lVert\widehat V^{-1}\rVertであり、 ∥V^eΛ^V^−1−eA∥≤∥E∥e∥A∥+∥E∥\lVert\widehat Ve^{\widehat\Lambda}\widehat V^{-1}-e^A\rVert\le\lVert E\rVert e^{\lVert A\rVert+\lVert E\rVert} が成り立つ。A≠0A\ne0ならば∥R∥∥V^−1∥/∥A∥=κ(V^)∥R∥/(∥A∥∥V^∥)\lVert R\rVert\lVert\widehat V^{-1}\rVert/\lVert A\rVert=\kappa(\widehat V)\lVert R\rVert/(\lVert A\rVert\lVert\widehat V\rVert)である。
  3. Rn\R^nのノルムが Euclid ノルムであり、V^\widehat Vが直交行列であるならば、∥E∥2=∥R∥2\lVert E\rVert_2=\lVert R\rVert_2かつκ2(V^)=1\kappa_2(\widehat V)=1である。AAが実対称行列ならば、直交行列QQと実対角行列Λ\Lambdaが存在して、V^=Q\widehat V=Q、Λ^=Λ\widehat\Lambda=\Lambdaに対してR=0R=0である。

証明.(1)を示す。A+E=A−(AV^−V^Λ^)V^−1=V^Λ^V^−1A+E=A-(A\widehat V-\widehat V\widehat\Lambda)\widehat V^{-1}=\widehat V\widehat\Lambda\widehat V^{-1}であるから、§E10.10 命題 4.1によりeA+E=V^eΛ^V^−1e^{A+E}=\widehat Ve^{\widehat\Lambda}\widehat V^{-1}である。各k∈N≥0k\in\NについてΛ^k=diag⁡(μ1k,…,μnk)\widehat\Lambda^k=\operatorname{diag}(\mu_1^k,\dots,\mu_n^k)であるから、行列指数関数の級数は成分ごとにeΛ^=diag⁡(eμ1,…,eμn)e^{\widehat\Lambda}=\operatorname{diag}(e^{\mu_1},\dots,e^{\mu_n})を与える。A=V^Λ^V^−1A=\widehat V\widehat\Lambda\widehat V^{-1}ならばR=0R=0であり、E=0E=0である。

(2)を示す。補題 1.1により∥E∥≤∥R∥∥V^−1∥\lVert E\rVert\le\lVert R\rVert\lVert\widehat V^{-1}\rVertである。(1)によりV^eΛ^V^−1−eA=eA+E−eA\widehat Ve^{\widehat\Lambda}\widehat V^{-1}-e^A=e^{A+E}-e^Aであり、補題 1.2をX=AX=Aに適用して二つ目の不等式を得る。V^\widehat Vは正則であるから∥V^∥>0\lVert\widehat V\rVert>0であり、κ(V^)=∥V^∥∥V^−1∥\kappa(\widehat V)=\lVert\widehat V\rVert\lVert\widehat V^{-1}\rVertから最後の等式が成り立つ。

(3)を示す。直交行列WWとh∈Rnh\in\R^nについて∥Wh∥2=∥h∥2\lVert Wh\rVert_2=\lVert h\rVert_2である。V^−1=V^T\widehat V^{-1}=\widehat V^{\mathsf T}は直交行列であるから、h↦V^Thh\mapsto\widehat V^{\mathsf T}hは Euclid ノルムの単位球面をそれ自身の上へ全単射に写し、∥E∥2=∥RV^T∥2=∥R∥2\lVert E\rVert_2=\lVert R\widehat V^{\mathsf T}\rVert_2=\lVert R\rVert_2である。同じ理由で∥V^∥2=∥V^−1∥2=1\lVert\widehat V\rVert_2=\lVert\widehat V^{-1}\rVert_2=1であり、κ2(V^)=1\kappa_2(\widehat V)=1である。AAが実対称行列ならば、§D3.15 定理 3.1により直交行列QQと実対角行列Λ\LambdaがQTAQ=ΛQ^{\mathsf T}AQ=\Lambdaを満たし、AQ−QΛ=Q(QTAQ−Λ)=0AQ-Q\Lambda=Q(Q^{\mathsf T}AQ-\Lambda)=0である。▨

例 1.4.R2\R^2にノルム∥⋅∥∞\lVert\cdot\rVert_\inftyを入れる。0<ε≤10<\varepsilon\le1、δ>0\delta>0とし、

A:=(010ε),V:=(110ε),Λ:=diag⁡(0,ε),Λ^:=diag⁡(0,ε+δ)A:=\begin{pmatrix}0&1\\0&\varepsilon\end{pmatrix},\qquad V:=\begin{pmatrix}1&1\\0&\varepsilon\end{pmatrix},\qquad\Lambda:=\operatorname{diag}(0,\varepsilon),\qquad\widehat\Lambda:=\operatorname{diag}(0,\varepsilon+\delta)

と置く。AV=VΛAV=V\Lambdaであり、命題 1.3 (1)により

eA=VeΛV−1=(1(eε−1)/ε0eε),V−1=(1−1/ε01/ε)e^A=Ve^\Lambda V^{-1}=\begin{pmatrix}1&(e^\varepsilon-1)/\varepsilon\\0&e^\varepsilon\end{pmatrix},\qquad V^{-1}=\begin{pmatrix}1&-1/\varepsilon\\0&1/\varepsilon\end{pmatrix}

である。固有ベクトルを正確に、固有値ε\varepsilonをε+δ\varepsilon+\deltaとして計算した結果を、V^:=V\widehat V:=VとΛ^\widehat\Lambdaで表す。R=AV−VΛ^=V(Λ−Λ^)R=AV-V\widehat\Lambda=V(\Lambda-\widehat\Lambda)の第1列は00、第2列は−δ(1,ε)T-\delta(1,\varepsilon)^{\mathsf T}であり、

E=−RV−1=(0δ/ε0δ)E=-RV^{-1}=\begin{pmatrix}0&\delta/\varepsilon\\0&\delta\end{pmatrix}

である。§E20.5 補題 3.5 (1)により∥R∥∞=δ\lVert R\rVert_\infty=\delta、∥V−1∥∞=1+1/ε\lVert V^{-1}\rVert_\infty=1+1/\varepsilon、ε≤1\varepsilon\le1から∥E∥∞=δ/ε\lVert E\rVert_\infty=\delta/\varepsilonであり、∥E∥∞<δ(1+1/ε)=∥R∥∞∥V−1∥∞\lVert E\rVert_\infty<\delta(1+1/\varepsilon)=\lVert R\rVert_\infty\lVert V^{-1}\rVert_\inftyである。∥V∥∞=2\lVert V\rVert_\infty=2であるからκ∞(V)=2(1+1/ε)\kappa_\infty(V)=2(1+1/\varepsilon)である。計算値とeAe^Aの差は

VeΛ^V−1−eA=(0eε(eδ−1)/ε0eε(eδ−1))Ve^{\widehat\Lambda}V^{-1}-e^A=\begin{pmatrix}0&e^\varepsilon(e^\delta-1)/\varepsilon\\0&e^\varepsilon(e^\delta-1)\end{pmatrix}

である。ε≤1\varepsilon\le1からeε−1≤(eε−1)/εe^\varepsilon-1\le(e^\varepsilon-1)/\varepsilonであり、∥VeΛ^V−1−eA∥∞=eε(eδ−1)/ε\lVert Ve^{\widehat\Lambda}V^{-1}-e^A\rVert_\infty=e^\varepsilon(e^\delta-1)/\varepsilon、∥eA∥∞=1+(eε−1)/ε\lVert e^A\rVert_\infty=1+(e^\varepsilon-1)/\varepsilonである。ε=10−8\varepsilon=10^{-8}、δ=10−16\delta=10^{-16}では、差のノルムは1.00000001…×10−81.00000001\ldots\times10^{-8}、∥eA∥∞=2.000000005…\lVert e^A\rVert_\infty=2.000000005\ldotsであり、相対誤差は5.00000003…×10−95.00000003\ldots\times10^{-9}である。他方∥A∥∞=1\lVert A\rVert_\infty=1であり、補題 1.2により、∥E′∥∞≤δ∥A∥∞\lVert E'\rVert_\infty\le\delta\lVert A\rVert_\inftyを満たす任意のE′∈M2(R)E'\in M_2(\R)に対して∥eA+E′−eA∥∞≤δe1+δ=2.718…×10−16\lVert e^{A+E'}-e^A\rVert_\infty\le\delta e^{1+\delta}=2.718\ldots\times10^{-16}である。対角化による計算値の誤差1.00000001…×10−81.00000001\ldots\times10^{-8}は、この上界の3.6×1073.6\times10^7倍より大きい。

2 Taylor 近似の打切り誤差

定理 2.1.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。m∈N≥0m\in\N、X∈Mn(R)X\in M_n(\R)とし、r:=∥X∥r:=\lVert X\rVert、Tm(X):=∑j=0mXj/j!T_m(X):=\sum_{j=0}^mX^j/j!と置く。このとき

∥eX−Tm(X)∥≤∑j=m+1∞rjj!≤errm+1(m+1)!\lVert e^X-T_m(X)\rVert\le\sum_{j=m+1}^\infty\frac{r^j}{j!}\le\frac{e^rr^{m+1}}{(m+1)!}

が成り立つ。

証明.§E10.10 命題 1.2により行列指数関数の級数は収束するので、eX−Tm(X)=lim⁡N→∞∑j=m+1NXj/j!e^X-T_m(X)=\lim_{N\to\infty}\sum_{j=m+1}^NX^j/j!である。補題 1.1により∥Xj∥≤rj\lVert X^j\rVert\le r^jであるから、各N>mN>mについて∥∑j=m+1NXj/j!∥≤∑j=m+1Nrj/j!\bigl\lVert\sum_{j=m+1}^NX^j/j!\bigr\rVert\le\sum_{j=m+1}^Nr^j/j!であり、ノルムの連続性により一つ目の不等式を得る。i∈N≥0i\in\Nについて(m+1+i)!/((m+1)! i!)(m+1+i)!/\bigl((m+1)!\,i!\bigr)は二項係数であって11以上であるから、

∑j=m+1∞rjj!=rm+1∑i=0∞ri(m+1+i)!≤rm+1(m+1)!∑i=0∞rii!=errm+1(m+1)!\sum_{j=m+1}^\infty\frac{r^j}{j!}=r^{m+1}\sum_{i=0}^\infty\frac{r^i}{(m+1+i)!}\le\frac{r^{m+1}}{(m+1)!}\sum_{i=0}^\infty\frac{r^i}{i!}=\frac{e^rr^{m+1}}{(m+1)!}

である。▨

3 引数の縮小と繰り返す二乗

定義 3.1.FFを浮動小数点数系、fl⁡\operatorname{fl}をFFの最近接丸めとし、n∈N≥1n\in\NN、X,Y∈Mn(F)X,Y\in M_n(F)とする。1≤i,k≤n1\le i,k\le nの各組について、fl⁡(xi1y1k),…,fl⁡(xinynk)\operatorname{fl}(x_{i1}y_{1k}),\dots,\operatorname{fl}(x_{in}y_{nk})を計算してこの順に逐次和で加える計算式を考える。n2n^2個の計算式をすべてFFとfl⁡\operatorname{fl}の浮動小数点算術で実行して得た値を(i,k)(i,k)成分とする行列を、XXとYYの積の 浮動小数点算術による計算値 (computed matrix product) といい、fl⁡(XY)\operatorname{fl}(XY)と書く。fl⁡(XY)\operatorname{fl}(XY)は、これらの計算式の実行が範囲条件を満たすときに考える。

補題 3.2.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。X,X^∈Mn(R)X,\widehat X\in M_n(\R)とδ≥0\delta\ge0が∥X^−X∥≤δ\lVert\widehat X-X\rVert\le\deltaを満たすとする。

  1. ∥X^2−X2∥≤2∥X∥δ+δ2\lVert\widehat X^2-X^2\rVert\le2\lVert X\rVert\delta+\delta^2が成り立つ。
  2. Rn\R^nのノルムを∥⋅∥∞\lVert\cdot\rVert_\inftyとする。FFを浮動小数点数系、uuをその単位丸め誤差、fl⁡\operatorname{fl}をFFの最近接丸めとし、nu<1nu<1とし、γn:=nu/(1−nu)\gamma_n:=nu/(1-nu)と置く。X^∈Mn(F)\widehat X\in M_n(F)であり、fl⁡(X^X^)\operatorname{fl}(\widehat X\widehat X)の計算が範囲条件を満たすならば ∥fl⁡(X^X^)−X^2∥∞≤γn∥X^∥∞2,∥fl⁡(X^X^)−X2∥∞≤2∥X∥∞δ+δ2+γn(∥X∥∞+δ)2\lVert\operatorname{fl}(\widehat X\widehat X)-\widehat X^2\rVert_\infty\le\gamma_n\lVert\widehat X\rVert_\infty^2,\qquad\lVert\operatorname{fl}(\widehat X\widehat X)-X^2\rVert_\infty\le2\lVert X\rVert_\infty\delta+\delta^2+\gamma_n(\lVert X\rVert_\infty+\delta)^2 が成り立つ。
  3. b≥∥X∥b\ge\lVert X\rVertとε≥0\varepsilon\ge0がδ≤bε\delta\le b\varepsilonを満たすとする。Z^:=X^2\widehat Z:=\widehat X^2かつγ:=0\gamma:=0とするか、(2)の仮定の下でZ^:=fl⁡(X^X^)\widehat Z:=\operatorname{fl}(\widehat X\widehat X)かつγ:=γn\gamma:=\gamma_nとする。このとき ∥Z^−X2∥≤b2((1+ε)2(1+γ)−1)\lVert\widehat Z-X^2\rVert\le b^2\bigl((1+\varepsilon)^2(1+\gamma)-1\bigr) が成り立つ。

証明.(1)を示す。E:=X^−XE:=\widehat X-Xと置くとX^2−X2=XE+EX+E2\widehat X^2-X^2=XE+EX+E^2であり、補題 1.1により∥X^2−X2∥≤2∥X∥∥E∥+∥E∥2≤2∥X∥δ+δ2\lVert\widehat X^2-X^2\rVert\le2\lVert X\rVert\lVert E\rVert+\lVert E\rVert^2\le2\lVert X\rVert\delta+\delta^2である。

(2)を示す。1≤i,k≤n1\le i,k\le nを固定する。fl⁡(X^X^)\operatorname{fl}(\widehat X\widehat X)の(i,k)(i,k)成分は、(x^i1,…,x^in)(\hat x_{i1},\dots,\hat x_{in})と(x^1k,…,x^nk)(\hat x_{1k},\dots,\hat x_{nk})の成分ごとの積を逐次和で加えた計算値であり、範囲条件により§E20.3 命題 5.1の仮定が満たされる。nu<1nu<1であるから、§E20.3 命題 5.1により

∣fl⁡(X^X^)ik−(X^2)ik∣≤γn∑j=1n∣x^ij∣∣x^jk∣=γn(∣X^∣∣X^∣)ik\bigl\lvert\operatorname{fl}(\widehat X\widehat X)_{ik}-(\widehat X^2)_{ik}\bigr\rvert\le\gamma_n\sum_{j=1}^n\lvert\hat x_{ij}\rvert\lvert\hat x_{jk}\rvert=\gamma_n\bigl(\lvert\widehat X\rvert\lvert\widehat X\rvert\bigr)_{ik}

である。ここで∣X^∣\lvert\widehat X\rvertはX^\widehat Xの成分の絶対値を並べた行列である。§E20.5 補題 3.5 (2)、補題 1.1、§E20.5 補題 3.5 (1)をこの順に用いると

∥fl⁡(X^X^)−X^2∥∞≤γn∥∣X^∣∣X^∣∥∞≤γn∥∣X^∣∥∞2=γn∥X^∥∞2\lVert\operatorname{fl}(\widehat X\widehat X)-\widehat X^2\rVert_\infty\le\gamma_n\bigl\lVert\lvert\widehat X\rvert\lvert\widehat X\rvert\bigr\rVert_\infty\le\gamma_n\bigl\lVert\lvert\widehat X\rvert\bigr\rVert_\infty^2=\gamma_n\lVert\widehat X\rVert_\infty^2

である。∥X^∥∞≤∥X∥∞+δ\lVert\widehat X\rVert_\infty\le\lVert X\rVert_\infty+\deltaであるから、この不等式と(1)と三角不等式から二つ目の不等式を得る。

(3)を示す。(1)または(2)により∥Z^−X2∥≤2∥X∥δ+δ2+γ(∥X∥+δ)2\lVert\widehat Z-X^2\rVert\le2\lVert X\rVert\delta+\delta^2+\gamma(\lVert X\rVert+\delta)^2である。右辺は∥X∥\lVert X\rVertとδ\deltaのそれぞれについて単調非減少であるから、∥X∥≤b\lVert X\rVert\le bとδ≤bε\delta\le b\varepsilonにより

∥Z^−X2∥≤b2(2ε+ε2+γ(1+ε)2)=b2((1+ε)2(1+γ)−1)\lVert\widehat Z-X^2\rVert\le b^2\bigl(2\varepsilon+\varepsilon^2+\gamma(1+\varepsilon)^2\bigr)=b^2\bigl((1+\varepsilon)^2(1+\gamma)-1\bigr)

である。▨

補題 3.3.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。X∈Mn(R)X\in M_n(\R)、s∈N≥0s\in\Nとし、r:=∥X∥r:=\lVert X\rVertと置く。Y^0∈Mn(R)\widehat Y_0\in M_n(\R)とδ0≥0\delta_0\ge0が∥Y^0−eX∥≤δ0\lVert\widehat Y_0-e^X\rVert\le\delta_0を満たすとし、ε0:=e−rδ0\varepsilon_0:=e^{-r}\delta_0と置く。

  1. Y^j+1:=Y^j2\widehat Y_{j+1}:=\widehat Y_j^2（0≤j<s0\le j<s）と置くと、0≤j≤s0\le j\le sを満たす各jjに対して ∥Y^j−e2jX∥≤e2jr((1+ε0)2j−1)\lVert\widehat Y_j-e^{2^jX}\rVert\le e^{2^jr}\bigl((1+\varepsilon_0)^{2^j}-1\bigr) が成り立つ。
  2. Rn\R^nのノルムを∥⋅∥∞\lVert\cdot\rVert_\inftyとし、FF、uu、fl⁡\operatorname{fl}、γn\gamma_nを補題 3.2 (2)のとおりとする。Y^0∈Mn(F)\widehat Y_0\in M_n(F)とし、Y^j+1:=fl⁡(Y^jY^j)\widehat Y_{j+1}:=\operatorname{fl}(\widehat Y_j\widehat Y_j)（0≤j<s0\le j<s）の各計算が範囲条件を満たすとする。このとき0≤j≤s0\le j\le sを満たす各jjに対して ∥Y^j−e2jX∥∞≤e2jr((1+ε0)2j(1+γn)2j−1−1)\lVert\widehat Y_j-e^{2^jX}\rVert_\infty\le e^{2^jr}\bigl((1+\varepsilon_0)^{2^j}(1+\gamma_n)^{2^j-1}-1\bigr) が成り立つ。

証明.(1)の場合にγ:=0\gamma:=0、(2)の場合にγ:=γn\gamma:=\gamma_nと置き、0≤j≤s0\le j\le sに対して

bj:=e2jr,εj:=(1+ε0)2j(1+γ)2j−1−1b_j:=e^{2^jr},\qquad\varepsilon_j:=(1+\varepsilon_0)^{2^j}(1+\gamma)^{2^j-1}-1

と置く。bj+1=bj2b_{j+1}=b_j^2と1+εj+1=(1+εj)2(1+γ)1+\varepsilon_{j+1}=(1+\varepsilon_j)^2(1+\gamma)が成り立ち、εj≥0\varepsilon_j\ge0である。b0ε0=δ0b_0\varepsilon_0=\delta_0であるから、仮定により∥Y^0−eX∥≤b0ε0\lVert\widehat Y_0-e^X\rVert\le b_0\varepsilon_0である。j<sj<sとし、∥Y^j−e2jX∥≤bjεj\lVert\widehat Y_j-e^{2^jX}\rVert\le b_j\varepsilon_jと仮定する。§E10.10 定理 2.2によりe2j+1X=e2jXe2jXe^{2^{j+1}X}=e^{2^jX}e^{2^jX}であり、補題 1.1により∥e2jX∥≤e∥2jX∥=bj\lVert e^{2^jX}\rVert\le e^{\lVert2^jX\rVert}=b_jである。補題 3.2 (3)を、XXをe2jXe^{2^jX}、X^\widehat XをY^j\widehat Y_j、δ\deltaをbjεjb_j\varepsilon_j、bbをbjb_j、ε\varepsilonをεj\varepsilon_jとして適用すると

∥Y^j+1−e2j+1X∥≤bj2((1+εj)2(1+γ)−1)=bj+1εj+1\lVert\widehat Y_{j+1}-e^{2^{j+1}X}\rVert\le b_j^2\bigl((1+\varepsilon_j)^2(1+\gamma)-1\bigr)=b_{j+1}\varepsilon_{j+1}

である。jjに関する帰納法により、0≤j≤s0\le j\le sを満たす各jjについて∥Y^j−e2jX∥≤bjεj\lVert\widehat Y_j-e^{2^jX}\rVert\le b_j\varepsilon_jが成り立ち、これは主張の不等式である。▨

定義 3.4.n∈N≥1n\in\NN、A∈Mn(R)A\in M_n(\R)、m,s∈N≥0m,s\in\Nとし、X:=2−sAX:=2^{-s}A、Tm(X):=∑j=0mXj/j!T_m(X):=\sum_{j=0}^mX^j/j!と置く。Y0:=Tm(X)Y_0:=T_m(X)、Yj+1:=Yj2Y_{j+1}:=Y_j^2（0≤j<s0\le j<s）で定まるYs=Tm(X)2sY_s=T_m(X)^{2^s}を、次数mmの Taylor 近似と縮小回数ssによるeAe^Aの 引数の縮小と繰り返す二乗 (scaling and squaring) の値という。FFを浮動小数点数系、fl⁡\operatorname{fl}をFFの最近接丸めとし、Y^0∈Mn(F)\widehat Y_0\in M_n(F)をTm(X)T_m(X)の近似とするとき、Y^j+1:=fl⁡(Y^jY^j)\widehat Y_{j+1}:=\operatorname{fl}(\widehat Y_j\widehat Y_j)（0≤j<s0\le j<s）で定まるY^s\widehat Y_sを、Y^0\widehat Y_0からの浮動小数点算術による値という。

定理 3.5.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。A∈Mn(R)A\in M_n(\R)、m,s∈N≥0m,s\in\Nとし、X:=2−sAX:=2^{-s}A、r:=∥X∥=2−s∥A∥r:=\lVert X\rVert=2^{-s}\lVert A\rVertと置く。

  1. 次が成り立つ。 ∥Tm(X)2s−eA∥≤e∥A∥((1+rm+1(m+1)!)2s−1)≤e∥A∥(exp⁡(∥A∥m+12sm(m+1)!)−1)\bigl\lVert T_m(X)^{2^s}-e^A\bigr\rVert\le e^{\lVert A\rVert}\Bigl(\Bigl(1+\frac{r^{m+1}}{(m+1)!}\Bigr)^{2^s}-1\Bigr)\le e^{\lVert A\rVert}\Bigl(\exp\Bigl(\frac{\lVert A\rVert^{m+1}}{2^{sm}(m+1)!}\Bigr)-1\Bigr)
  2. Rn\R^nのノルムを∥⋅∥∞\lVert\cdot\rVert_\inftyとし、FF、uu、fl⁡\operatorname{fl}、γn\gamma_nを補題 3.2 (2)のとおりとする。Y^0∈Mn(F)\widehat Y_0\in M_n(F)とη≥0\eta\ge0が∥Y^0−Tm(X)∥∞≤η\lVert\widehat Y_0-T_m(X)\rVert_\infty\le\etaを満たし、定義 3.4のY^1,…,Y^s\widehat Y_1,\dots,\widehat Y_sの各計算が範囲条件を満たすとする。ε0:=rm+1/(m+1)!+e−rη\varepsilon_0:=r^{m+1}/(m+1)!+e^{-r}\etaと置くと ∥Y^s−eA∥∞≤e∥A∥∞((1+ε0)2s(1+γn)2s−1−1)\lVert\widehat Y_s-e^A\rVert_\infty\le e^{\lVert A\rVert_\infty}\bigl((1+\varepsilon_0)^{2^s}(1+\gamma_n)^{2^s-1}-1\bigr) が成り立つ。

証明.定理 2.1により∥Tm(X)−eX∥≤errm+1/(m+1)!\lVert T_m(X)-e^X\rVert\le e^rr^{m+1}/(m+1)!である。2sX=A2^sX=Aであるからe2sX=eAe^{2^sX}=e^Aかつe2sr=e∥A∥e^{2^sr}=e^{\lVert A\rVert}である。補題 3.3 (1)をY^0:=Tm(X)\widehat Y_0:=T_m(X)、δ0:=errm+1/(m+1)!\delta_0:=e^rr^{m+1}/(m+1)!、j:=sj:=sとして適用すると、Y^s=Tm(X)2s\widehat Y_s=T_m(X)^{2^s}であり、ε0=rm+1/(m+1)!\varepsilon_0=r^{m+1}/(m+1)!であるから、一つ目の不等式を得る。x≥0x\ge0とN∈N≥0N\in\Nについて(1+x)N≤eNx(1+x)^N\le e^{Nx}であり、2srm+1=2−sm∥A∥m+12^sr^{m+1}=2^{-sm}\lVert A\rVert^{m+1}であるから、二つ目の不等式が成り立つ。(2)の仮定の下では、三角不等式により∥Y^0−eX∥∞≤errm+1/(m+1)!+η=erε0\lVert\widehat Y_0-e^X\rVert_\infty\le e^rr^{m+1}/(m+1)!+\eta=e^r\varepsilon_0である。補題 3.3 (2)をδ0:=erε0\delta_0:=e^r\varepsilon_0、j:=sj:=sとして適用すると、(2)の不等式を得る。▨

4 有理近似

補題 4.1.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。M∈Mn(R)M\in M_n(\R)が∥M∥<1\lVert M\rVert<1を満たすならば、I−MI-Mは正則であり、∥(I−M)−1∥≤(1−∥M∥)−1\lVert(I-M)^{-1}\rVert\le(1-\lVert M\rVert)^{-1}が成り立つ。

証明.h∈Rnh\in\R^nをとる。§E20.2 補題 1.2 (1)により∥h∥≤∥(I−M)h∥+∥Mh∥≤∥(I−M)h∥+∥M∥∥h∥\lVert h\rVert\le\lVert(I-M)h\rVert+\lVert Mh\rVert\le\lVert(I-M)h\rVert+\lVert M\rVert\lVert h\rVertであるから、(1−∥M∥)∥h∥≤∥(I−M)h∥(1-\lVert M\rVert)\lVert h\rVert\le\lVert(I-M)h\rVertである。(I−M)h=0(I-M)h=0ならばh=0h=0であるから、正方行列I−MI-Mは単射であり、正則である。任意のk∈Rnk\in\R^nについて、h:=(I−M)−1kh:=(I-M)^{-1}kに上の不等式を適用すると∥(I−M)−1k∥≤(1−∥M∥)−1∥k∥\lVert(I-M)^{-1}k\rVert\le(1-\lVert M\rVert)^{-1}\lVert k\rVertである。▨

定義 4.2.n∈N≥1n\in\NN、X∈Mn(R)X\in M_n(\R)とし、p,qp,qを実係数の多項式とする。q(X)q(X)が正則であるとき、q(X)Z=p(X)q(X)Z=p(X)を満たすただ一つの行列Z=q(X)−1p(X)∈Mn(R)Z=q(X)^{-1}p(X)\in M_n(\R)を、組(p,q)(p,q)によるXXの 有理関数値 (rational function of a matrix) という。組(1+x/2, 1−x/2)(1+x/2,\,1-x/2)は指数関数の[1/1][1/1]型の Padé 近似を与える組であり（§E20.17 命題 2.1 (1)）、I−X/2I-X/2が正則であるとき、この組によるXXの有理関数値をR1,1(X):=(I−X/2)−1(I+X/2)R_{1,1}(X):=(I-X/2)^{-1}(I+X/2)と書く。

命題 4.3.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。X∈Mn(R)X\in M_n(\R)とし、r:=∥X∥<2r:=\lVert X\rVert<2とする。

  1. I−X/2I-X/2は正則であり、∥(I−X/2)−1∥≤(1−r/2)−1\lVert(I-X/2)^{-1}\rVert\le(1-r/2)^{-1}が成り立つ。特にR1,1(X)R_{1,1}(X)が定まる。
  2. 次が成り立つ。 ∥eX−R1,1(X)∥≤r3er12(1−r/2)\lVert e^X-R_{1,1}(X)\rVert\le\frac{r^3e^r}{12(1-r/2)}

証明.補題 4.1をM:=X/2M:=X/2に適用すると、∥M∥=r/2<1\lVert M\rVert=r/2<1であるから(1)が成り立つ。

(2)を示す。G:=(I−X/2)eX−(I+X/2)G:=(I-X/2)e^X-(I+X/2)と置く。行列指数関数の級数に左からI−X/2I-X/2を掛けてXXの冪ごとにまとめると、(I−X/2)eX(I-X/2)e^XのXkX^kの係数は、k=0k=0で11、k≥1k\ge1で1/k!−1/(2(k−1)!)=(2−k)/(2⋅k!)1/k!-1/\bigl(2(k-1)!\bigr)=(2-k)/(2\cdot k!)であり、k=1,2k=1,2ではそれぞれ1/21/2、00である。したがって

G=∑k=3∞2−k2⋅k!XkG=\sum_{k=3}^\infty\frac{2-k}{2\cdot k!}X^k

である。k≥3k\ge3ならばk(k−1)≥6k(k-1)\ge6であるから

k−22⋅k!=12k(k−1)(k−3)!≤112(k−3)!\frac{k-2}{2\cdot k!}=\frac1{2k(k-1)(k-3)!}\le\frac1{12(k-3)!}

であり、補題 1.1により

∥G∥≤∑k=3∞rk12(k−3)!=r3er12\lVert G\rVert\le\sum_{k=3}^\infty\frac{r^k}{12(k-3)!}=\frac{r^3e^r}{12}

である。eX−R1,1(X)=(I−X/2)−1Ge^X-R_{1,1}(X)=(I-X/2)^{-1}Gであるから、(1)と補題 1.1により主張の不等式が成り立つ。▨

注意 4.4.定義 4.2の記号で、q(X)q(X)が正則であるとする。Z=q(X)−1p(X)Z=q(X)^{-1}p(X)の第kk列zkz_kは、p(X)p(X)の第kk列ckc_kを右辺とするq(X)zk=ckq(X)z_k=c_kのただ一つの解である。§E20.5 定理 2.5 (1)により、q(X)q(X)の部分ピボット付き消去は厳密算術でどの段でも失敗せず、ピボット付き LU 分解(P,L,U)(P,L,U)を与える。§E20.5 定理 2.5 (2)により、各kkについてLy=PckLy=Pc_kの前進代入とUz=yUz=yの後退代入がzkz_kを与える。p(X)p(X)とq(X)q(X)を得た後のこの計算はq(X)−1q(X)^{-1}を作らず、その四則演算の回数は、§E20.5 定理 2.5 (3)により分解の2n3/3−n2/2−n/62n^3/3-n^2/2-n/6回とnn組の代入のn(2n2−n)n(2n^2-n)回の和8n3/3−3n2/2−n/68n^3/3-3n^2/2-n/6である。

系 4.5.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。A∈Mn(R)A\in M_n(\R)、s∈N≥0s\in\Nとし、X:=2−sAX:=2^{-s}A、r:=2−s∥A∥r:=2^{-s}\lVert A\rVertと置いて、r<2r<2とする。このとき

∥R1,1(X)2s−eA∥≤e∥A∥((1+r312(1−r/2))2s−1)\bigl\lVert R_{1,1}(X)^{2^s}-e^A\bigr\rVert\le e^{\lVert A\rVert}\Bigl(\Bigl(1+\frac{r^3}{12(1-r/2)}\Bigr)^{2^s}-1\Bigr)

が成り立つ。

証明.命題 4.3 (2)により∥R1,1(X)−eX∥≤err3/(12(1−r/2))\lVert R_{1,1}(X)-e^X\rVert\le e^rr^3/\bigl(12(1-r/2)\bigr)である。補題 3.3 (1)をY^0:=R1,1(X)\widehat Y_0:=R_{1,1}(X)、δ0:=err3/(12(1−r/2))\delta_0:=e^rr^3/\bigl(12(1-r/2)\bigr)、j:=sj:=sとして適用し、e2sX=eAe^{2^sX}=e^Aとe2sr=e∥A∥e^{2^sr}=e^{\lVert A\rVert}を用いると、主張の不等式を得る。▨

5 過渡的な増大と導関数

例 5.1.K>1K>1とし、

A:=(−1K0−1),P:=diag⁡(K,1),J:=(−110−1)A:=\begin{pmatrix}-1&K\\0&-1\end{pmatrix},\qquad P:=\operatorname{diag}(K,1),\qquad J:=\begin{pmatrix}-1&1\\0&-1\end{pmatrix}

と置く。AAの固有値は−1-1だけであり、A=PJP−1A=PJP^{-1}である。§E10.10 定理 5.1により、t∈Rt\in\Rに対して

etA=PetJP−1=e−tP(1t01)P−1=e−t(1Kt01)e^{tA}=Pe^{tJ}P^{-1}=e^{-t}P\begin{pmatrix}1&t\\0&1\end{pmatrix}P^{-1}=e^{-t}\begin{pmatrix}1&Kt\\0&1\end{pmatrix}

である。R2\R^2にノルム∥⋅∥∞\lVert\cdot\rVert_\inftyを入れると、§E20.5 補題 3.5 (1)により、t≥0t\ge0ならば∥etA∥∞=e−t(1+Kt)\lVert e^{tA}\rVert_\infty=e^{-t}(1+Kt)である。ddt(e−t(1+Kt))=e−t(K−1−Kt)\frac{d}{dt}\bigl(e^{-t}(1+Kt)\bigr)=e^{-t}(K-1-Kt)は0≤t<1−1/K0\le t<1-1/Kで正、t>1−1/Kt>1-1/Kで負であるから、∥etA∥∞\lVert e^{tA}\rVert_\inftyはt≥0t\ge0においてt=1−1/Kt=1-1/Kで最大値Ke1/K−1Ke^{1/K-1}をとる。K=100K=100では最大値は37.157…37.157\ldotsであり、固有値−1-1から定まるe−t≤1e^{-t}\le1を上回る。y0:=(0,1)Ty_0:=(0,1)^{\mathsf T}に対するy′=Ayy'=Ay、y(0)=y0y(0)=y_0の解は、§E10.10 定理 3.1によりy(t)=e−t(Kt,1)Ty(t)=e^{-t}(Kt,1)^{\mathsf T}であり、その第1成分はt=1t=1で最大値K/eK/eをとる。定理 3.5の上界の因子e∥A∥∞=eK+1e^{\lVert A\rVert_\infty}=e^{K+1}は、∥eA∥∞=(1+K)/e\lVert e^A\rVert_\infty=(1+K)/eのeK+2/(1+K)e^{K+2}/(1+K)倍である。

命題 5.2.n∈N≥1n\in\NN、λ∈R\lambda\in\Rとし、N∈Mn(R)N\in M_n(\R)がN≠0N\ne0とN2=0N^2=0を満たすとして、A:=λI+NA:=\lambda I+Nと置く。

  1. 任意の実係数の多項式ppに対してp(A)=p(λ)I+p′(λ)Np(A)=p(\lambda)I+p'(\lambda)Nが成り立つ。またeA=eλ(I+N)e^A=e^\lambda(I+N)が成り立つ。
  2. AAの最小多項式は(t−λ)2(t-\lambda)^2である。ffをλ\lambdaを含む開区間で定義された実数値関数でλ\lambdaで微分可能なものとし、f(A)f(A)を§E3.30 命題 3.1により定めると、f(A)=f(λ)I+f′(λ)Nf(A)=f(\lambda)I+f'(\lambda)Nが成り立つ。

証明.(1)を示す。ppの次数をddとすると、多項式の Taylor 展開により、変数xxの多項式としてp(λ+x)=∑k=0dp(k)(λ)xk/k!p(\lambda+x)=\sum_{k=0}^dp^{(k)}(\lambda)x^k/k!である。x↦Nx\mapsto Nで定まる代入はλ+x\lambda+xをAAへ写す環準同型であるから、p(A)=∑k=0dp(k)(λ)Nk/k!p(A)=\sum_{k=0}^dp^{(k)}(\lambda)N^k/k!であり、k≥2k\ge2でNk=0N^k=0であるからp(A)=p(λ)I+p′(λ)Np(A)=p(\lambda)I+p'(\lambda)Nである。λI\lambda IとNNは可換であるから、§E10.10 命題 2.1によりeA=eλIeNe^A=e^{\lambda I}e^Nである。級数の定義によりeλI=eλIe^{\lambda I}=e^\lambda Iであり、N2=0N^2=0からeN=I+Ne^N=I+Nである。

(2)を示す。(A−λI)2=N2=0(A-\lambda I)^2=N^2=0であるから、AAの最小多項式は(t−λ)2(t-\lambda)^2を割り切る。A−λI=N≠0A-\lambda I=N\ne0であるから最小多項式はt−λt-\lambdaでなく、n≥1n\ge1であるから定数でもないので、(t−λ)2(t-\lambda)^2である。p(t):=f(λ)+f′(λ)(t−λ)p(t):=f(\lambda)+f'(\lambda)(t-\lambda)はp(λ)=f(λ)p(\lambda)=f(\lambda)とp′(λ)=f′(λ)p'(\lambda)=f'(\lambda)を満たすので、§E3.30 命題 3.1によりf(A)=p(A)f(A)=p(A)であり、(1)によりf(A)=f(λ)I+f′(λ)Nf(A)=f(\lambda)I+f'(\lambda)Nである。▨

例 5.3.R2\R^2にノルム∥⋅∥∞\lVert\cdot\rVert_\inftyを入れる。K>0K>0とし、A:=N:=(0K00)A:=N:=\begin{pmatrix}0&K\\0&0\end{pmatrix}と置く。AAの固有値は00だけであり、命題 5.2 (1)によりeA=I+Ne^A=I+N、∥eA∥∞=1+K\lVert e^A\rVert_\infty=1+Kである。

  1. p(x):=1+(e−1)xp(x):=1+(e-1)xはx=0x=0とx=1x=1でexe^xに一致し、固有値00でp(0)=e0p(0)=e^0を満たす。命題 5.2 (1)によりp(A)−eA=(p′(0)−1)N=(e−2)Np(A)-e^A=(p'(0)-1)N=(e-2)Nであり、∥p(A)−eA∥∞=(e−2)K=0.71828…K\lVert p(A)-e^A\rVert_\infty=(e-2)K=0.71828\ldots Kである。
  2. h>0h>0とし、fh(x):=ex+hsin⁡(x/h)f_h(x):=e^x+h\sin(x/h)と置く。任意の実数xxに対して∣fh(x)−ex∣≤h\lvert f_h(x)-e^x\rvert\le hである。fh′(0)=2f_h'(0)=2であるから、命題 5.2 (2)と命題 5.2 (1)によりfh(A)−eA=(fh(0)−1)I+(fh′(0)−1)N=Nf_h(A)-e^A=(f_h(0)-1)I+(f_h'(0)-1)N=Nであり、∥fh(A)−eA∥∞=K\lVert f_h(A)-e^A\rVert_\infty=Kはhhによらない。sup⁡x∈R∣fh(x)−ex∣≤h\sup_{x\in\R}\lvert f_h(x)-e^x\rvert\le hであるが、∥fh(A)−eA∥∞\lVert f_h(A)-e^A\rVert_\inftyはh→+0h\to+0で00に収束しない。

6 線形微分方程式系への適用

系 6.1.n∈N≥1n\in\NNとし、Rn\R^nにノルムを固定して、Mn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書く。A∈Mn(R)A\in M_n(\R)、y0∈Rn∖{0}y_0\in\R^n\setminus\{0\}、t∈Rt\in\Rとし、yyをy′=Ayy'=Ay、y(0)=y0y(0)=y_0の解とする。

  1. Y^∈Mn(R)\widehat Y\in M_n(\R)とΔ≥0\Delta\ge0が∥Y^−etA∥≤Δ\lVert\widehat Y-e^{tA}\rVert\le\Deltaを満たすとし、y^:=Y^y0\hat y:=\widehat Yy_0と置く。このときy(t)≠0y(t)\ne0であり、 ∥y^−y(t)∥≤Δ∥y0∥,∥y^−y(t)∥∥y(t)∥≤Δ∥e−tA∥≤Δe∣t∣∥A∥\lVert\hat y-y(t)\rVert\le\Delta\lVert y_0\rVert,\qquad\frac{\lVert\hat y-y(t)\rVert}{\lVert y(t)\rVert}\le\Delta\lVert e^{-tA}\rVert\le\Delta e^{\lvert t\rvert\lVert A\rVert} が成り立つ。
  2. m,s∈N≥0m,s\in\Nとし、y^:=Tm(2−stA)2sy0\hat y:=T_m(2^{-s}tA)^{2^s}y_0と置くと ∥y^−y(t)∥∥y(t)∥≤e2∣t∣∥A∥(exp⁡((∣t∣∥A∥)m+12sm(m+1)!)−1)\frac{\lVert\hat y-y(t)\rVert}{\lVert y(t)\rVert}\le e^{2\lvert t\rvert\lVert A\rVert}\Bigl(\exp\Bigl(\frac{(\lvert t\rvert\lVert A\rVert)^{m+1}}{2^{sm}(m+1)!}\Bigr)-1\Bigr) が成り立つ。

証明.§E10.10 定理 3.1をt0=0t_0=0、η=y0\boldsymbol\eta=y_0に適用するとy(t)=etAy0y(t)=e^{tA}y_0である。§E10.10 定理 2.2によりetAe^{tA}は正則で(etA)−1=e−tA(e^{tA})^{-1}=e^{-tA}であるから、y(t)≠0y(t)\ne0かつy0=e−tAy(t)y_0=e^{-tA}y(t)である。y^−y(t)=(Y^−etA)y0\hat y-y(t)=(\widehat Y-e^{tA})y_0であるから、§E20.2 補題 1.2 (1)により∥y^−y(t)∥≤Δ∥y0∥≤Δ∥e−tA∥∥y(t)∥\lVert\hat y-y(t)\rVert\le\Delta\lVert y_0\rVert\le\Delta\lVert e^{-tA}\rVert\lVert y(t)\rVertであり、補題 1.1により∥e−tA∥≤e∣t∣∥A∥\lVert e^{-tA}\rVert\le e^{\lvert t\rvert\lVert A\rVert}である。これで(1)は示された。定理 3.5 (1)をtAtAに適用すると、Y^:=Tm(2−stA)2s\widehat Y:=T_m(2^{-s}tA)^{2^s}とΔ:=e∣t∣∥A∥(exp⁡((∣t∣∥A∥)m+1/(2sm(m+1)!))−1)\Delta:=e^{\lvert t\rvert\lVert A\rVert}\bigl(\exp\bigl((\lvert t\rvert\lVert A\rVert)^{m+1}/(2^{sm}(m+1)!)\bigr)-1\bigr)は(1)の仮定を満たし、その二つ目の不等式から(2)を得る。▨

例 6.2.R2\R^2にノルム∥⋅∥∞\lVert\cdot\rVert_\inftyを入れ、例 5.1のAAをK=10K=10として、y0:=(0,1)Ty_0:=(0,1)^{\mathsf T}、t=1t=1、m=8m=8とする。基準解は例 5.1の公式による

eA=e−1(11001),y(1)=e−1(101)=(3.6787944117…0.36787944117…)e^A=e^{-1}\begin{pmatrix}1&10\\0&1\end{pmatrix},\qquad y(1)=e^{-1}\begin{pmatrix}10\\1\end{pmatrix}=\begin{pmatrix}3.6787944117\ldots\\0.36787944117\ldots\end{pmatrix}

であり、∥A∥∞=11\lVert A\rVert_\infty=11、∥eA∥∞=11/e=4.0466…\lVert e^A\rVert_\infty=11/e=4.0466\ldots、∥e−A∥∞=11e=29.901…\lVert e^{-A}\rVert_\infty=11e=29.901\ldotsである。次の表の (a) は定理 3.5 (1)の一つ目の上界、(b) はT8(2−sA)2sT_8(2^{-s}A)^{2^s}を有理数の厳密な演算で計算した値の誤差∥T8(2−sA)2s−eA∥∞\lVert T_8(2^{-s}A)^{2^s}-e^A\rVert_\inftyである。(c) は定理 3.5 (2)の上界であり、F=F(2,53,−1022,1023)F=F(2,53,-1022,1023)、u=2−53u=2^{-53}、最近接偶数丸めとし、Y^0\widehat Y_0をT8(2−sA)T_8(2^{-s}A)の各成分を丸めた行列、η:=u∥T8(2−sA)∥∞\eta:=u\lVert T_8(2^{-s}A)\rVert_\inftyとする。T8(2−sA)T_8(2^{-s}A)の成分は00であるか正規範囲にあるので、§E20.1 定理 2.2 (2)により∥Y^0−T8(2−sA)∥∞≤η\lVert\widehat Y_0-T_8(2^{-s}A)\rVert_\infty\le\etaである。(d) は CPython の float でY^s\widehat Y_sを定義 3.1の順序で計算した値の誤差∥Y^s−eA∥∞\lVert\widehat Y_s-e^A\rVert_\inftyであり、観察である。表の各ssについて、二乗の各計算は範囲条件を満たす。上界と誤差は 80 桁の十進演算で評価した。

ss r=11/2sr=11/2^s (a) (b) (c) (d)
00 1111 3.8905×1083.8905\times10^{8} 2.2549×10−42.2549\times10^{-4} 3.8905×1083.8905\times10^{8} 2.2549×10−42.2549\times10^{-4}
22 2.752.75 6.1609×1036.1609\times10^{3} 1.6132×10−91.6132\times10^{-9} 6.1609×1036.1609\times10^{3} 1.6132×10−91.6132\times10^{-9}
44 0.68750.6875 9.0584×10−29.0584\times10^{-2} 2.0366×10−142.0366\times10^{-14} 9.0584×10−29.0584\times10^{-2} 2.0150×10−142.0150\times10^{-14}
66 0.171880.17188 1.3822×10−61.3822\times10^{-6} 2.9638×10−192.9638\times10^{-19} 1.3834×10−61.3834\times10^{-6} 8.9075×10−158.9075\times10^{-15}
88 0.0429690.042969 2.1091×10−112.1091\times10^{-11} 4.4691×10−244.4691\times10^{-24} 5.0985×10−95.0985\times10^{-9} 8.9075×10−158.9075\times10^{-15}
1010 0.0107420.010742 3.2182×10−163.2182\times10^{-16} 6.7992×10−296.7992\times10^{-29} 2.0394×10−82.0394\times10^{-8} 1.2470×10−131.2470\times10^{-13}
1212 0.00268550.0026855 4.9106×10−214.9106\times10^{-21} 1.0367×10−331.0367\times10^{-33} 8.1656×10−88.1656\times10^{-8} 5.6729×10−135.6729\times10^{-13}

各行で (b) は (a) 以下、(d) は (c) 以下である。表の範囲で (a) と (b) はssについて減少し、(c) はs=8s=8で最小である。(d) はs=10,12s=10,12でs=8s=8の値より大きい。(a) と (c) は因子e∥A∥∞=e11=59874.1…e^{\lVert A\rVert_\infty}=e^{11}=59874.1\ldotsを含み、この因子は∥eA∥∞\lVert e^A\rVert_\inftyの1.4×1041.4\times10^4倍より大きい。s=8s=8のy^:=Y^8y0\hat y:=\widehat Y_8y_0はY^8\widehat Y_8の第2列であり、この積に丸めは入らない。その相対誤差∥y^−y(1)∥∞/∥y(1)∥∞\lVert\hat y-y(1)\rVert_\infty/\lVert y(1)\rVert_\inftyは2.2067×10−152.2067\times10^{-15}（観察）であり、系 6.1 (1)を (c) の値Δ=5.0985×10−9\Delta=5.0985\times10^{-9}に適用した上界はΔ∥e−A∥∞=1.5245×10−7\Delta\lVert e^{-A}\rVert_\infty=1.5245\times10^{-7}である。

7 演習

問題 7.1.補題 1.2の証明を完成させよ。

解答.

k∈N≥1k\in\NNとする。0≤j≤k−10\le j\le k-1について(X+E)j+1Xk−1−j−(X+E)jXk−j=(X+E)jEXk−1−j(X+E)^{j+1}X^{k-1-j}-(X+E)^jX^{k-j}=(X+E)^jEX^{k-1-j}であり、これをj=0,…,k−1j=0,\dots,k-1について加えると

(X+E)k−Xk=∑j=0k−1(X+E)jEXk−1−j(X+E)^k-X^k=\sum_{j=0}^{k-1}(X+E)^jEX^{k-1-j}

である。補題 1.1により、右辺の各項のノルムは∥X+E∥j∥E∥∥X∥k−1−j≤∥E∥(∥X∥+∥E∥)k−1\lVert X+E\rVert^j\lVert E\rVert\lVert X\rVert^{k-1-j}\le\lVert E\rVert(\lVert X\rVert+\lVert E\rVert)^{k-1}以下であるから、∥(X+E)k−Xk∥≤k∥E∥(∥X∥+∥E∥)k−1\lVert(X+E)^k-X^k\rVert\le k\lVert E\rVert(\lVert X\rVert+\lVert E\rVert)^{k-1}である。§E10.10 命題 1.2によりeX+Ee^{X+E}とeXe^Xの級数はどちらも収束するので、eX+E−eX=∑k=1∞((X+E)k−Xk)/k!e^{X+E}-e^X=\sum_{k=1}^\infty\bigl((X+E)^k-X^k\bigr)/k!であり、ノルムの連続性により

∥eX+E−eX∥≤∥E∥∑k=1∞(∥X∥+∥E∥)k−1(k−1)!=∥E∥e∥X∥+∥E∥\lVert e^{X+E}-e^X\rVert\le\lVert E\rVert\sum_{k=1}^\infty\frac{(\lVert X\rVert+\lVert E\rVert)^{k-1}}{(k-1)!}=\lVert E\rVert e^{\lVert X\rVert+\lVert E\rVert}

である。▨

問題 7.2.n∈N≥1n\in\NN、λ∈R∖{2}\lambda\in\R\setminus\{2\}とし、N∈Mn(R)N\in M_n(\R)がN≠0N\ne0とN2=0N^2=0を満たすとして、A:=λI+NA:=\lambda I+Nと置く。r1,1(x):=(2+x)/(2−x)r_{1,1}(x):=(2+x)/(2-x)とする。I−A/2I-A/2が正則であり、R1,1(A)=r1,1(λ)I+r1,1′(λ)NR_{1,1}(A)=r_{1,1}(\lambda)I+r_{1,1}'(\lambda)Nが成り立つことを示せ。さらにλ=1\lambda=1のときR1,1(A)−eAR_{1,1}(A)-e^Aを求めよ。また、Rn\R^nにノルムを固定してMn(R)M_n(\R)の元の作用素ノルムを∥⋅∥\lVert\cdot\rVertと書くとき、∥N∥≥2+∣λ∣\lVert N\rVert\ge2+\lvert\lambda\rvertならば∥A∥≥2\lVert A\rVert\ge2であることを示せ。

解答.

c:=1−λ/2c:=1-\lambda/2と置くとλ≠2\lambda\ne2からc≠0c\ne0であり、I−A/2=cI−N/2I-A/2=cI-N/2である。a:=c−1a:=c^{-1}、b:=c−2/2b:=c^{-2}/2と置くと、N2=0N^2=0により

(cI−N/2)(aI+bN)=acI+(bc−a/2)N=I(cI-N/2)(aI+bN)=acI+(bc-a/2)N=I

であるから、I−A/2I-A/2は正則であり、その逆行列はaI+bNaI+bNである。I+A/2=(1+λ/2)I+N/2I+A/2=(1+\lambda/2)I+N/2であるから、N2=0N^2=0により

R1,1(A)=(aI+bN)((1+λ/2)I+N/2)=a(1+λ/2)I+(a/2+b(1+λ/2))NR_{1,1}(A)=(aI+bN)\bigl((1+\lambda/2)I+N/2\bigr)=a(1+\lambda/2)I+\bigl(a/2+b(1+\lambda/2)\bigr)N

である。a(1+λ/2)=(2+λ)/(2−λ)=r1,1(λ)a(1+\lambda/2)=(2+\lambda)/(2-\lambda)=r_{1,1}(\lambda)であり、

a2+b(1+λ2)=c−22((1−λ2)+(1+λ2))=c−2=4(2−λ)2\frac a2+b\Bigl(1+\frac\lambda2\Bigr)=\frac{c^{-2}}2\Bigl(\bigl(1-\frac\lambda2\bigr)+\bigl(1+\frac\lambda2\bigr)\Bigr)=c^{-2}=\frac4{(2-\lambda)^2}

である。r1,1′(x)=((2−x)+(2+x))/(2−x)2=4/(2−x)2r_{1,1}'(x)=\bigl((2-x)+(2+x)\bigr)/(2-x)^2=4/(2-x)^2であるから、R1,1(A)=r1,1(λ)I+r1,1′(λ)NR_{1,1}(A)=r_{1,1}(\lambda)I+r_{1,1}'(\lambda)Nが成り立つ。λ=1\lambda=1ではr1,1(1)=3r_{1,1}(1)=3、r1,1′(1)=4r_{1,1}'(1)=4であり、命題 5.2 (1)によりeA=e(I+N)e^A=e(I+N)であるから

R1,1(A)−eA=(3−e)I+(4−e)N=0.28171…I+1.28171…NR_{1,1}(A)-e^A=(3-e)I+(4-e)N=0.28171\ldots I+1.28171\ldots N

である。Rn\R^nのノルムに関する作用素ノルムについて、∥N∥=∥A−λI∥≤∥A∥+∣λ∣\lVert N\rVert=\lVert A-\lambda I\rVert\le\lVert A\rVert+\lvert\lambda\rvertであるから、∥N∥≥2+∣λ∣\lVert N\rVert\ge2+\lvert\lambda\rvertならば∥A∥≥2\lVert A\rVert\ge2である。▨

前提記事

9 本の記事・単元を表示