§E20.31保存構造を考慮した時間積分

最終更新

Rd\R^dの開集合上のC2C^2級関数KK、VVから定まる Hamilton 系q′=∇K(p)q'=\nabla K(p)、p′=−∇V(q)p'=-\nabla V(q)の解に沿って、H(q,p):=K(p)+V(q)H(q,p):=K(p)+V(q)の値は一定である。調和振動子d=1d=1、K(p)=p2/2K(p)=p^2/2、V(q)=q2/2V(q)=q^2/2を刻み1/101/10の前進 Euler 法で時刻100100まで解くと、(1,0)(1,0)から出た近似値のHHの値は初期値の2×1042\times10^4倍を超える。

∇K(p)\nabla K(p)と∇V(q)\nabla V(q)の一方だけを残した系の解は陽に書くことができ、これを部分流という。部分流は symplectic 写像であり、symplectic 写像の合成もまた symplectic 写像であるから、部分流を刻みh/2h/2、hh、h/2h/2で交互に合成した Störmer–Verlet 法の一歩は symplectic 写像である。前進 Euler 法と同じ調和振動子、刻み、初期値では、Störmer–Verlet 法の近似値のHHの値はすべての歩で初期値での値との差が0.001250.00125以下にとどまる。この評価は、一歩がHHと異なる二次形式(1−h2/4)q2+p2(1-h^2/4)q^2+p^2を正確に保つことから従う。

本記事は、部分流の合成として symplectic Euler 法と Störmer–Verlet 法を構成し、調和振動子と振り子でそれらの振る舞いを調べる。

1 部分流と symplectic 写像

定義 1.1.d∈N≥1d\in\NNとし、IdI_dをdd次単位行列として

J:=(0Id−Id0)∈R2d×2dJ:=\begin{pmatrix}0&I_d\\-I_d&0\end{pmatrix}\in\R^{2d\times2d}

と置く。開集合W⊆R2dW\subseteq\R^{2d}上のC1C^1級写像Φ ⁣:W→R2d\Phi\colon W\to\R^{2d}が、任意のz∈Wz\in Wについて

DΦ(z)TJDΦ(z)=JD\Phi(z)^{\mathsf T}JD\Phi(z)=J

を満たすとき、Φ\Phiを symplectic 写像 (symplectic map) という。ここでDΦ(z)∈R2d×2dD\Phi(z)\in\R^{2d\times2d}はΦ\Phiのzzにおける Jacobi 行列である。

定義 1.2.d∈N≥1d\in\NNとし、Uq,Up⊆RdU_q,U_p\subseteq\R^dを開集合、U:=Uq×Up⊆R2dU:=U_q\times U_p\subseteq\R^{2d}とする。R2d\R^{2d}の元をz=(q,p)z=(q,p)（q,p∈Rdq,p\in\R^d）と書く。K∈C2(Up;R)K\in C^2(U_p;\R)とV∈C2(Uq;R)V\in C^2(U_q;\R)に対してH(q,p):=K(p)+V(q)H(q,p):=K(p)+V(q)と置き、HHを 分離形 Hamilton 関数 (separable Hamiltonian)、UUに値をとるC1C^1級写像z=(q,p)z=(q,p)に対する連立常微分方程式

q′=∇K(p),p′=−∇V(q)q'=\nabla K(p),\qquad p'=-\nabla V(q)

をHHの Hamilton 系 (Hamiltonian system) という。h∈Rh\in\Rに対して

WhV:={(q,p)∈U∣p−h∇V(q)∈Up},WhK:={(q,p)∈U∣q+h∇K(p)∈Uq}W^V_h:=\{(q,p)\in U\mid p-h\nabla V(q)\in U_p\},\qquad W^K_h:=\{(q,p)\in U\mid q+h\nabla K(p)\in U_q\}

と置き、写像

φhV ⁣:WhV→U, (q,p)↦(q, p−h∇V(q)),φhK ⁣:WhK→U, (q,p)↦(q+h∇K(p), p)\varphi^V_h\colon W^V_h\to U,\ (q,p)\mapsto\bigl(q,\ p-h\nabla V(q)\bigr),\qquad \varphi^K_h\colon W^K_h\to U,\ (q,p)\mapsto\bigl(q+h\nabla K(p),\ p\bigr)

を定める。φhV\varphi^V_hを刻みhhの キック (kick)、φhK\varphi^K_hを刻みhhの ドリフト (drift) といい、両者をあわせてHHの 部分流 (split flows) という。

補題 1.3.定義 1.2の設定で、z0=(q0,p0)∈Uz_0=(q_0,p_0)\in Uとし、I⊆RI\subseteq\Rを00を含む区間とする。任意のs∈Is\in Iについて(q0,p0−s∇V(q0))∈U(q_0,p_0-s\nabla V(q_0))\in Uであるならば、写像I→UI\to U、s↦φsV(z0)s\mapsto\varphi^V_s(z_0)はK=0K=0とした Hamilton 系q′=0q'=0、p′=−∇V(q)p'=-\nabla V(q)の解である。任意のs∈Is\in Iについて(q0+s∇K(p0),p0)∈U(q_0+s\nabla K(p_0),p_0)\in Uであるならば、写像I→UI\to U、s↦φsK(z0)s\mapsto\varphi^K_s(z_0)はV=0V=0とした Hamilton 系q′=∇K(p)q'=\nabla K(p)、p′=0p'=0の解である。

JJを定義 1.1の行列とする。任意の(q,p)∈U(q,p)\in Uについて∇H(q,p)=(∇V(q),∇K(p))\nabla H(q,p)=(\nabla V(q),\nabla K(p))であり、区間上で定義されUUに値をとるC1C^1級写像zzが Hamilton 系を満たすこととz′=J∇H(z)z'=J\nabla H(z)を満たすことは同値である。I0⊆RI_0\subseteq\Rを区間とし、C1C^1級写像z=(q,p) ⁣:I0→Uz=(q,p)\colon I_0\to Uが Hamilton 系を満たすならば、H∘zH\circ zはI0I_0上で一定である。

証明.s∈Is\in IについてφsV(z0)=(q0,p0−s∇V(q0))\varphi^V_s(z_0)=(q_0,p_0-s\nabla V(q_0))はssの一次式であるから、写像s↦φsV(z0)s\mapsto\varphi^V_s(z_0)はII上でC1C^1級であり、その導関数は(0,−∇V(q0))(0,-\nabla V(q_0))である。φsV(z0)\varphi^V_s(z_0)の第一成分はq0q_0であるから、この導関数はK=0K=0とした Hamilton 系の右辺(0,−∇V(q))(0,-\nabla V(q))に(q,p)=φsV(z0)(q,p)=\varphi^V_s(z_0)を代入した値に等しい。したがって写像s↦φsV(z0)s\mapsto\varphi^V_s(z_0)はK=0K=0とした Hamilton 系の解である。s∈Is\in IについてφsK(z0)=(q0+s∇K(p0),p0)\varphi^K_s(z_0)=(q_0+s\nabla K(p_0),p_0)はssの一次式であるから、写像s↦φsK(z0)s\mapsto\varphi^K_s(z_0)はII上でC1C^1級であり、その導関数は(∇K(p0),0)(\nabla K(p_0),0)である。φsK(z0)\varphi^K_s(z_0)の第二成分はp0p_0であるから、この導関数はV=0V=0とした Hamilton 系の右辺(∇K(p),0)(\nabla K(p),0)に(q,p)=φsK(z0)(q,p)=\varphi^K_s(z_0)を代入した値に等しい。したがって写像s↦φsK(z0)s\mapsto\varphi^K_s(z_0)はV=0V=0とした Hamilton 系の解である。

(q,p)∈U(q,p)\in Uとする。H(q,p)=K(p)+V(q)H(q,p)=K(p)+V(q)のqqに関する勾配は∇V(q)\nabla V(q)、ppに関する勾配は∇K(p)\nabla K(p)であるから、∇H(q,p)=(∇V(q),∇K(p))\nabla H(q,p)=(\nabla V(q),\nabla K(p))である。x,y∈Rdx,y\in\R^dについてJ(x,y)=(y,−x)J(x,y)=(y,-x)であるから、J∇H(q,p)=(∇K(p),−∇V(q))J\nabla H(q,p)=(\nabla K(p),-\nabla V(q))である。したがって、区間上で定義されUUに値をとるC1C^1級写像z=(q,p)z=(q,p)について、z′=J∇H(z)z'=J\nabla H(z)を満たすことはq′=∇K(p)q'=\nabla K(p)とp′=−∇V(q)p'=-\nabla V(q)を満たすことと同値である。

I0I_0の内部をI0∘I_0^\circと書き、τ∈I0∘\tau\in I_0^\circとする。§E4.3 定理 1.1をq∣I0∘q|_{I_0^\circ}とVV、p∣I0∘p|_{I_0^\circ}とKKに適用すると

(H∘z)′(τ)=∇V(q(τ))⋅q′(τ)+∇K(p(τ))⋅p′(τ)=∇V(q(τ))⋅∇K(p(τ))−∇K(p(τ))⋅∇V(q(τ))=0(H\circ z)'(\tau)=\nabla V(q(\tau))\cdot q'(\tau)+\nabla K(p(\tau))\cdot p'(\tau)=\nabla V(q(\tau))\cdot\nabla K(p(\tau))-\nabla K(p(\tau))\cdot\nabla V(q(\tau))=0

である。a,b∈I0a,b\in I_0、a<ba<bとすると、[a,b]⊆I0[a,b]\subseteq I_0かつ(a,b)⊆I0∘(a,b)\subseteq I_0^\circであり、H∘zH\circ zは[a,b][a,b]上で連続で(a,b)(a,b)上で導関数00をもつ。平均値の定理によりH(z(a))=H(z(b))H(z(a))=H(z(b))である。▨

補題 1.4.定義 1.2の設定で、次が成り立つ。

  1. 任意のh∈Rh\in\Rについて、WhVW^V_hとWhKW^K_hは開集合であり、φhV\varphi^V_hとφhK\varphi^K_hは symplectic 写像である。
  2. 任意のh∈Rh\in\RについてφhV(WhV)=W−hV\varphi^V_h(W^V_h)=W^V_{-h}であり、任意のz∈WhVz\in W^V_hについてφ−hV(φhV(z))=z\varphi^V_{-h}(\varphi^V_h(z))=zである。φhK\varphi^K_hについても同じことが成り立つ。
  3. m∈N≥1m\in\NN、σ1,…,σm∈{K,V}\sigma_1,\dots,\sigma_m\in\{K,V\}、h1,…,hm∈Rh_1,\dots,h_m\in\Rとする。z∈Uz\in Uに対してz0:=zz_0:=zと置き、zj−1∈Whjσjz_{j-1}\in W^{\sigma_j}_{h_j}である限りzj:=φhjσj(zj−1)z_j:=\varphi^{\sigma_j}_{h_j}(z_{j-1})と定める。z1,…,zmz_1,\dots,z_mがすべて定まるzzの全体をWWとし、S(z):=zmS(z):=z_mと置く。列(σm,…,σ1)(\sigma_m,\dots,\sigma_1)と(−hm,…,−h1)(-h_m,\dots,-h_1)から同じ方法で定まる集合と写像をW−W^-、S−S^-とする。このときWWは開集合、S ⁣:W→R2dS\colon W\to\R^{2d}は symplectic 写像であり、S(W)=W−S(W)=W^-かつ任意のz∈Wz\in WについてS−(S(z))=zS^-(S(z))=zが成り立つ。

証明.(1)を示す。h∈Rh\in\Rとする。V∈C2V\in C^2であるから∇V ⁣:Uq→Rd\nabla V\colon U_q\to\R^dはC1C^1級であり、写像(q,p)↦p−h∇V(q)(q,p)\mapsto p-h\nabla V(q)はUU上で連続である。WhVW^V_hはこの写像による開集合UpU_pの逆像であるから開集合であり、φhV\varphi^V_hはC1C^1級である。(q,p)∈WhV(q,p)\in W^V_hとし、A:=D(∇V)(q)A:=D(\nabla V)(q)をVVの Hesse 行列とすると、§E4.4 定理 2.1によりAT=AA^{\mathsf T}=Aである。M:=DφhV(q,p)M:=D\varphi^V_h(q,p)と定義 1.1のJJについて

M=(Id0−hAId),MTJM=(Id−hAT0Id)(−hAId−Id0)=(h(AT−A)Id−Id0)=JM=\begin{pmatrix}I_d&0\\-hA&I_d\end{pmatrix},\qquad M^{\mathsf T}JM=\begin{pmatrix}I_d&-hA^{\mathsf T}\\0&I_d\end{pmatrix}\begin{pmatrix}-hA&I_d\\-I_d&0\end{pmatrix}=\begin{pmatrix}h(A^{\mathsf T}-A)&I_d\\-I_d&0\end{pmatrix}=J

である。同様にWhKW^K_hは開集合、φhK\varphi^K_hはC1C^1級であり、(q,p)∈WhK(q,p)\in W^K_hについてB:=D(∇K)(p)B:=D(\nabla K)(p)はBT=BB^{\mathsf T}=Bを満たし、N:=DφhK(q,p)N:=D\varphi^K_h(q,p)について

N=(IdhB0Id),NTJN=(Id0hBTId)(0Id−Id−hB)=(0Id−Idh(BT−B))=JN=\begin{pmatrix}I_d&hB\\0&I_d\end{pmatrix},\qquad N^{\mathsf T}JN=\begin{pmatrix}I_d&0\\hB^{\mathsf T}&I_d\end{pmatrix}\begin{pmatrix}0&I_d\\-I_d&-hB\end{pmatrix}=\begin{pmatrix}0&I_d\\-I_d&h(B^{\mathsf T}-B)\end{pmatrix}=J

である。

(2)を示す。h∈Rh\in\R、z=(q,p)∈WhVz=(q,p)\in W^V_hとし、w:=φhV(z)=(q,p−h∇V(q))w:=\varphi^V_h(z)=(q,p-h\nabla V(q))と置く。w∈Uw\in Uであり、wwの第二成分にh∇V(q)h\nabla V(q)を加えたものはp∈Upp\in U_pであるから、w∈W−hVw\in W^V_{-h}かつφ−hV(w)=z\varphi^V_{-h}(w)=zである。したがってφhV(WhV)⊆W−hV\varphi^V_h(W^V_h)\subseteq W^V_{-h}である。hhを−h-hに替えると、任意のw∈W−hVw\in W^V_{-h}についてφ−hV(w)∈WhV\varphi^V_{-h}(w)\in W^V_hかつφhV(φ−hV(w))=w\varphi^V_h(\varphi^V_{-h}(w))=wであるから、W−hV⊆φhV(WhV)W^V_{-h}\subseteq\varphi^V_h(W^V_h)である。φhK\varphi^K_hについては、第一成分にh∇K(p)h\nabla K(p)を加える操作について同じ議論が成り立つ。

(3)を示す。m=1m=1のとき、主張は(1)と(2)である。m≥2m\ge2とし、長さm−1m-1の任意の列について主張が成り立つとする。最初のm−1m-1項の列(σ1,…,σm−1)(\sigma_1,\dots,\sigma_{m-1})と(h1,…,hm−1)(h_1,\dots,h_{m-1})から定まる集合と写像をW′W'、S′S'、その逆順の列から定まるものをW′−W'^-、S′−S'^-とする。Φ:=φhmσm\Phi:=\varphi^{\sigma_m}_{h_m}と置くと、

W={z∈W′∣S′(z)∈Whmσm},S=Φ∘S′∣WW=\{z\in W'\mid S'(z)\in W^{\sigma_m}_{h_m}\},\qquad S=\Phi\circ S'|_W

である。S′S'は連続でありWhmσmW^{\sigma_m}_{h_m}は開集合であるから、WWは開集合である。§E4.3 定理 1.1により、z∈Wz\in WについてDS(z)=DΦ(S′(z))DS′(z)DS(z)=D\Phi(S'(z))DS'(z)であり、右辺はzzについて連続であるからSSはC1C^1級である。さらに

DS(z)TJDS(z)=DS′(z)T(DΦ(S′(z))TJDΦ(S′(z)))DS′(z)=DS′(z)TJDS′(z)=JDS(z)^{\mathsf T}JDS(z)=DS'(z)^{\mathsf T}\bigl(D\Phi(S'(z))^{\mathsf T}JD\Phi(S'(z))\bigr)DS'(z)=DS'(z)^{\mathsf T}JDS'(z)=J

である。逆順の列は(σm,−hm)(\sigma_m,-h_m)から始まり、その後に最初のm−1m-1項の逆順の列が続くから、

W−={w∈W−hmσm∣φ−hmσm(w)∈W′−},S−(w)=S′−(φ−hmσm(w))W^-=\{w\in W^{\sigma_m}_{-h_m}\mid\varphi^{\sigma_m}_{-h_m}(w)\in W'^-\},\qquad S^-(w)=S'^-\bigl(\varphi^{\sigma_m}_{-h_m}(w)\bigr)

である。z∈Wz\in Wとし、w:=S(z)=Φ(S′(z))w:=S(z)=\Phi(S'(z))と置く。(2)によりw∈W−hmσmw\in W^{\sigma_m}_{-h_m}かつφ−hmσm(w)=S′(z)\varphi^{\sigma_m}_{-h_m}(w)=S'(z)であり、帰納法の仮定によりS′(z)∈W′−S'(z)\in W'^-かつS′−(S′(z))=zS'^-(S'(z))=zである。したがってw∈W−w\in W^-かつS−(w)=zS^-(w)=zである。以上により、長さmmの任意の列についてS(W)⊆W−S(W)\subseteq W^-とS−∘S=id⁡WS^-\circ S=\id_Wが成り立つ。これを逆順の列に適用すると、その逆順の列はもとの列であるから、S−(W−)⊆WS^-(W^-)\subseteq WとS∘S−=id⁡W−S\circ S^-=\id_{W^-}が成り立つ。w∈W−w\in W^-はw=S(S−(w))w=S(S^-(w))と書くことができるからW−⊆S(W)W^-\subseteq S(W)であり、S(W)=W−S(W)=W^-である。▨

2 symplectic Euler 法と Störmer–Verlet 法

定義 2.1.定義 1.2の設定でh∈Rh\in\Rとする。補題 1.4 (3)の構成を列(V,K)(V,K)と(h,h)(h,h)、列(K,V)(K,V)と(h,h)(h,h)、列(V,K,V)(V,K,V)と(h/2,h,h/2)(h/2,h,h/2)に適用して得られる写像を、それぞれ

ΦhVK:=φhK∘φhV,ΦhKV:=φhV∘φhK,Sh:=φh/2V∘φhK∘φh/2V\Phi^{VK}_h:=\varphi^K_h\circ\varphi^V_h,\qquad\Phi^{KV}_h:=\varphi^V_h\circ\varphi^K_h,\qquad S_h:=\varphi^V_{h/2}\circ\varphi^K_h\circ\varphi^V_{h/2}

とし、それぞれの定義域をWhVKW^{VK}_h、WhKVW^{KV}_h、WhSVW^{\mathrm{SV}}_hとする。(q1,p1)(q_1,p_1)をそれぞれの写像による(q,p)(q,p)の像とすると、

ΦhVK:p1=p−h∇V(q),q1=q+h∇K(p1),ΦhKV:q1=q+h∇K(p),p1=p−h∇V(q1),Sh:pˉ=p−h2∇V(q),q1=q+h∇K(pˉ),p1=pˉ−h2∇V(q1)\begin{aligned} \Phi^{VK}_h&:\quad p_1=p-h\nabla V(q),\quad q_1=q+h\nabla K(p_1),\\ \Phi^{KV}_h&:\quad q_1=q+h\nabla K(p),\quad p_1=p-h\nabla V(q_1),\\ S_h&:\quad\bar p=p-\tfrac h2\nabla V(q),\quad q_1=q+h\nabla K(\bar p),\quad p_1=\bar p-\tfrac h2\nabla V(q_1) \end{aligned}

である。Ω:=R×U\Omega:=\R\times U上の右辺f(t,q,p):=(∇K(p),−∇V(q))f(t,q,p):=(\nabla K(p),-\nabla V(q))に対して、D:={(t,z,h)∈Ω×(0,∞)∣z∈WhVK}D:=\{(t,z,h)\in\Omega\times(0,\infty)\mid z\in W^{VK}_h\}上の増分関数(t,z,h)↦(ΦhVK(z)−z)/h(t,z,h)\mapsto(\Phi^{VK}_h(z)-z)/hを対応させる一段法と、ΦhKV\Phi^{KV}_hから同じ方法で定まる一段法を symplectic Euler 法 (symplectic Euler method) といい、ShS_hから同じ方法で定まる一段法を Störmer–Verlet 法 (Störmer–Verlet method) という。h≤0h\le0を含む任意のh∈Rh\in\Rについて、ΦhVK\Phi^{VK}_h、ΦhKV\Phi^{KV}_h、ShS_hをそれぞれの一歩写像と呼ぶ。

命題 2.2.定義 1.2の設定でUp=RdU_p=\R^dとし、h∈Rh\in\Rとすると、次が成り立つ。

  1. 任意のa,b∈Ra,b\in\Rについて、φaV\varphi^V_aの定義域はUUであり、φaV∘φbV=φa+bV\varphi^V_a\circ\varphi^V_b=\varphi^V_{a+b}である。
  2. n∈N≥1n\in\NN、z∈Uz\in Uとし、w:=φ−h/2V(z)w:=\varphi^V_{-h/2}(z)と置く。z,Sh(z),…,Shn−1(z)z,S_h(z),\dots,S_h^{n-1}(z)がすべてWhSVW^{\mathrm{SV}}_hに属することと、w,ΦhVK(w),…,(ΦhVK)n−1(w)w,\Phi^{VK}_h(w),\dots,(\Phi^{VK}_h)^{n-1}(w)がすべてWhVKW^{VK}_hに属することは同値であり、そのときShn(z)=φh/2V((ΦhVK)n(w))S_h^n(z)=\varphi^V_{h/2}\bigl((\Phi^{VK}_h)^n(w)\bigr)である。

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

定義 2.3.m∈N≥1m\in\NNとし、各h∈Rh\in\Rに対して集合Wh⊆RmW_h\subseteq\R^mと単射Ψh ⁣:Wh→Rm\Psi_h\colon W_h\to\R^mが与えられているとする。Ψ−h(W−h)\Psi_{-h}(W_{-h})上で定まる写像Ψh∗:=(Ψ−h)−1\Psi^*_h:=(\Psi_{-h})^{-1}の族(Ψh∗)h∈R(\Psi^*_h)_{h\in\R}を(Ψh)h∈R(\Psi_h)_{h\in\R}の 随伴 (adjoint) という。任意のh∈Rh\in\RについてΨ−h(W−h)=Wh\Psi_{-h}(W_{-h})=W_hかつΨh∗=Ψh\Psi^*_h=\Psi_hが成り立つとき、(Ψh)h∈R(\Psi_h)_{h\in\R}は 時間対称 (time-symmetric) であるという。

補題 2.4.m,n∈N≥1m,n\in\NN、r∈N≥0r\in\Nとし、Rm\R^mとRn\R^nに Euclid ノルムを入れる。U⊆RmU\subseteq\R^mを開集合、g∈Cr+1(U;Rn)g\in C^{r+1}(U;\R^n)とし、a∈Ua\in Uとv∈Rmv\in\R^mが{a+sv∣0≤s≤1}⊆U\{a+sv\mid0\le s\le1\}\subseteq Uを満たすとする。M≥0M\ge0が存在して、0≤s≤10\le s\le1について∥Dr+1g(a+sv)∥op≤M\|D^{r+1}g(a+sv)\|_{\mathrm{op}}\le Mが成り立つならば、

∥g(a+v)−∑j=0r1j!Djg(a)[v,…,v]∥≤n M(r+1)!∥v∥r+1\Bigl\|g(a+v)-\sum_{j=0}^{r}\frac1{j!}D^jg(a)[v,\dots,v]\Bigr\|\le\frac{\sqrt n\,M}{(r+1)!}\|v\|^{r+1}

が成り立つ。m=1m=1のとき、Djg(a)[v,…,v]=vjg(j)(a)D^jg(a)[v,\dots,v]=v^jg^{(j)}(a)であり、∥Djg(a)∥op=∥g(j)(a)∥\|D^jg(a)\|_{\mathrm{op}}=\|g^{(j)}(a)\|である。

証明.1≤i≤n1\le i\le nについてeie_iをRn\R^nの第ii標準基底ベクトルとし、gi:=ei⋅gg_i:=e_i\cdot gと置く。gig_iはggと線形写像w↦ei⋅ww\mapsto e_i\cdot wの合成であるから、§E4.3 定理 1.1によりgi∈Cr+1(U;R)g_i\in C^{r+1}(U;\R)であり、0≤j≤r+10\le j\le r+1についてDjgi(x)[v1,…,vj]=ei⋅Djg(x)[v1,…,vj]D^jg_i(x)[v_1,\dots,v_j]=e_i\cdot D^jg(x)[v_1,\dots,v_j]である。∣ei⋅w∣≤∥w∥|e_i\cdot w|\le\|w\|であるから、0≤s≤10\le s\le1について∥Dr+1gi(a+sv)∥op≤M\|D^{r+1}g_i(a+sv)\|_{\mathrm{op}}\le Mである。§E4.5 定理 1.1と§E4.5 系 1.2をgig_iに適用すると、左辺のベクトルの第ii成分RiR_iは∣Ri∣≤M∥v∥r+1/(r+1)!|R_i|\le M\|v\|^{r+1}/(r+1)!を満たす。∥(R1,…,Rn)∥≤nmax⁡i∣Ri∣\|(R_1,\dots,R_n)\|\le\sqrt n\max_i|R_i|から主張を得る。

m=1m=1のとき、全微分Dg(a)Dg(a)はs↦sg′(a)s\mapsto sg'(a)であり、jjについての帰納法によりDjg(a)[s1,…,sj]=s1⋯sjg(j)(a)D^jg(a)[s_1,\dots,s_j]=s_1\cdots s_jg^{(j)}(a)である。∣s1∣,…,∣sj∣≤1|s_1|,\dots,|s_j|\le1にわたる上限をとると∥Djg(a)∥op=∥g(j)(a)∥\|D^jg(a)\|_{\mathrm{op}}=\|g^{(j)}(a)\|を得る。▨

命題 2.5.定義 1.2の設定で、Rd\R^dとR2d=Rd×Rd\R^{2d}=\R^d\times\R^dに Euclid ノルムを入れる。

  1. 任意のh∈Rh\in\Rについて、ΦhVK\Phi^{VK}_hとΦhKV\Phi^{KV}_hは symplectic 写像であり、Φ−hVK(W−hVK)=WhKV\Phi^{VK}_{-h}(W^{VK}_{-h})=W^{KV}_hかつ任意のz∈W−hVKz\in W^{VK}_{-h}についてΦhKV(Φ−hVK(z))=z\Phi^{KV}_h(\Phi^{VK}_{-h}(z))=zである。すなわち(ΦhVK)h∈R(\Phi^{VK}_h)_{h\in\R}の随伴は(ΦhKV)h∈R(\Phi^{KV}_h)_{h\in\R}である。
  2. d=1d=1、Uq=Up=RU_q=U_p=\R、K(p)=p2/2K(p)=p^2/2、V(q)=q2/2V(q)=q^2/2とする。任意のh≠0h\ne0についてΦ−hVK∘ΦhVK≠id⁡R2\Phi^{VK}_{-h}\circ\Phi^{VK}_h\ne\id_{\R^2}かつΦ−hKV∘ΦhKV≠id⁡R2\Phi^{KV}_{-h}\circ\Phi^{KV}_h\ne\id_{\R^2}であり、(ΦhVK)h∈R(\Phi^{VK}_h)_{h\in\R}と(ΦhKV)h∈R(\Phi^{KV}_h)_{h\in\R}はいずれも時間対称でない。
  3. t0<Tt_0<Tとし、[t0,T][t_0,T]を含む開区間JJと、Hamilton 系を満たすC1C^1級写像z=(q,p) ⁣:J→Uz=(q,p)\colon J\to Uを取る。ρ>0\rho>0とM>0M>0が次を満たすとする。任意のτ∈[t0,T]\tau\in[t_0,T]について閉球Bˉ(q(τ),ρ)\bar B(q(\tau),\rho)はUqU_qに、Bˉ(p(τ),ρ)\bar B(p(\tau),\rho)はUpU_pに含まれ、これらの閉球のτ∈[t0,T]\tau\in[t_0,T]にわたる和集合をそれぞれNqN_q、NpN_pとすると、任意のx∈Nqx\in N_qとy∈Npy\in N_pについて ∥∇V(x)∥≤M,∥D(∇V)(x)∥op≤M,∥∇K(y)∥≤M,∥D(∇K)(y)∥op≤M\|\nabla V(x)\|\le M,\quad\|D(\nabla V)(x)\|_{\mathrm{op}}\le M,\quad\|\nabla K(y)\|\le M,\quad\|D(\nabla K)(y)\|_{\mathrm{op}}\le M が成り立つ。γ:=d M\gamma:=\sqrt d\,Mと置く。このとき、t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{ρ/M,T−t}0<h\le\min\{\rho/M,T-t\}を満たす任意のt,ht,hについてz(t)∈WhVK∩WhKVz(t)\in W^{VK}_h\cap W^{KV}_hであり、Φh\Phi_hがΦhVK\Phi^{VK}_hとΦhKV\Phi^{KV}_hのいずれであっても ∥z(t+h)−Φh(z(t))∥≤2γMh2\|z(t+h)-\Phi_h(z(t))\|\le2\gamma Mh^2 が成り立つ。
  4. (2)のK,VK,Vについて、z(τ):=(cos⁡τ,−sin⁡τ)z(\tau):=(\cos\tau,-\sin\tau)とw(τ):=(sin⁡τ,cos⁡τ)w(\tau):=(\sin\tau,\cos\tau)はR\R上で Hamilton 系を満たし、任意のh>0h>0について ∥z(h)−ΦhVK(z(0))∥≥h22,∥w(h)−ΦhKV(w(0))∥≥h22\|z(h)-\Phi^{VK}_h(z(0))\|\ge\frac{h^2}2,\qquad\|w(h)-\Phi^{KV}_h(w(0))\|\ge\frac{h^2}2 が成り立つ。したがって、t0=0t_0=0、T=1T=1とするとき、どのようなC≥0C\ge0とh0>0h_0>0に対しても、0<h≤min⁡{h0,1}0<h\le\min\{h_0,1\}を満たすすべてのhhについて∥z(h)−ΦhVK(z(0))∥≤Ch3\|z(h)-\Phi^{VK}_h(z(0))\|\le Ch^3が成り立つことはなく、wwとΦhKV\Phi^{KV}_hについても同様である。

証明.(1)を示す。ΦhVK\Phi^{VK}_hとΦhKV\Phi^{KV}_hが symplectic 写像であることは補題 1.4 (3)による。補題 1.4 (3)を列(V,K)(V,K)と(−h,−h)(-h,-h)に適用すると、S=Φ−hVKS=\Phi^{VK}_{-h}、W=W−hVKW=W^{VK}_{-h}であり、逆順の列(K,V)(K,V)と(h,h)(h,h)から定まる写像と集合はS−=ΦhKVS^-=\Phi^{KV}_h、W−=WhKVW^-=W^{KV}_hである。したがってΦ−hVK(W−hVK)=WhKV\Phi^{VK}_{-h}(W^{VK}_{-h})=W^{KV}_hかつΦhKV∘Φ−hVK=id⁡W−hVK\Phi^{KV}_h\circ\Phi^{VK}_{-h}=\id_{W^{VK}_{-h}}であり、Φ−hVK\Phi^{VK}_{-h}は単射で、その逆写像はΦhKV\Phi^{KV}_hである。

(2)を示す。∇V(q)=q\nabla V(q)=q、∇K(p)=p\nabla K(p)=pであるから、すべての部分流の定義域はR2\R^2であり、z∈R2z\in\R^2を列ベクトルとみなすと

ΦhVK(z)=(1−h2h−h1)z,ΦhKV(z)=(1h−h1−h2)z\Phi^{VK}_h(z)=\begin{pmatrix}1-h^2&h\\-h&1\end{pmatrix}z,\qquad\Phi^{KV}_h(z)=\begin{pmatrix}1&h\\-h&1-h^2\end{pmatrix}z

である。Φ−hVK∘ΦhVK\Phi^{VK}_{-h}\circ\Phi^{VK}_hの行列の(2,2)(2,2)成分はh⋅h+1⋅1=1+h2h\cdot h+1\cdot1=1+h^2、Φ−hKV∘ΦhKV\Phi^{KV}_{-h}\circ\Phi^{KV}_hの行列の(1,1)(1,1)成分は1⋅1+(−h)(−h)=1+h21\cdot1+(-h)(-h)=1+h^2であり、h≠0h\ne0ならばいずれも11でない。族が時間対称ならば、Ψh∗=Ψh\Psi^*_h=\Psi_hによりWhW_h上でΨ−h∘Ψh=id⁡\Psi_{-h}\circ\Psi_h=\idとなるから、二つの族は時間対称でない。

(3)を示す。t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{ρ/M,T−t}0<h\le\min\{\rho/M,T-t\}を取る。∇K\nabla Kと∇V\nabla VはC1C^1級であるから、§E4.3 定理 1.1によりq′=∇K∘pq'=\nabla K\circ pとp′=−∇V∘qp'=-\nabla V\circ qはJJ上でC1C^1級であり、τ∈J\tau\in Jについて

q′′(τ)=−D(∇K)(p(τ))∇V(q(τ)),p′′(τ)=−D(∇V)(q(τ))∇K(p(τ))q''(\tau)=-D(\nabla K)(p(\tau))\nabla V(q(\tau)),\qquad p''(\tau)=-D(\nabla V)(q(\tau))\nabla K(p(\tau))

である。τ∈[t0,T]\tau\in[t_0,T]についてq(τ)∈Nqq(\tau)\in N_q、p(τ)∈Npp(\tau)\in N_pであるから∥q′′(τ)∥≤M2\|q''(\tau)\|\le M^2、∥p′′(τ)∥≤M2\|p''(\tau)\|\le M^2である。[t,t+h]⊆[t0,T][t,t+h]\subseteq[t_0,T]であるから、補題 2.4をm=1m=1、n=dn=d、r=1r=1でqqとppに適用すると

∥q(t+h)−q(t)−h∇K(p(t))∥≤γM2h2,∥p(t+h)−p(t)+h∇V(q(t))∥≤γM2h2\|q(t+h)-q(t)-h\nabla K(p(t))\|\le\frac{\gamma M}2h^2,\qquad\|p(t+h)-p(t)+h\nabla V(q(t))\|\le\frac{\gamma M}2h^2

である。Bˉ(q(t),ρ)\bar B(q(t),\rho)とBˉ(p(t),ρ)\bar B(p(t),\rho)は凸であり、それぞれNqN_q、NpN_pに含まれるから、補題 2.4をr=0r=0で∇V\nabla Vと∇K\nabla Kに適用すると、x,x′∈Bˉ(q(t),ρ)x,x'\in\bar B(q(t),\rho)とy,y′∈Bˉ(p(t),ρ)y,y'\in\bar B(p(t),\rho)について

∥∇V(x)−∇V(x′)∥≤γ∥x−x′∥,∥∇K(y)−∇K(y′)∥≤γ∥y−y′∥\|\nabla V(x)-\nabla V(x')\|\le\gamma\|x-x'\|,\qquad\|\nabla K(y)-\nabla K(y')\|\le\gamma\|y-y'\|

である。

Φh=ΦhVK\Phi_h=\Phi^{VK}_hの場合、p1:=p(t)−h∇V(q(t))p_1:=p(t)-h\nabla V(q(t))は∥p1−p(t)∥≤hM≤ρ\|p_1-p(t)\|\le hM\le\rhoを満たすからp1∈Bˉ(p(t),ρ)⊆Upp_1\in\bar B(p(t),\rho)\subseteq U_pであり、q1:=q(t)+h∇K(p1)q_1:=q(t)+h\nabla K(p_1)は∥∇K(p1)∥≤M\|\nabla K(p_1)\|\le Mから∥q1−q(t)∥≤hM≤ρ\|q_1-q(t)\|\le hM\le\rhoを満たすからq1∈Uqq_1\in U_qである。したがってz(t)∈WhVKz(t)\in W^{VK}_hかつΦhVK(z(t))=(q1,p1)\Phi^{VK}_h(z(t))=(q_1,p_1)である。∥∇K(p1)−∇K(p(t))∥≤γ∥p1−p(t)∥≤γMh\|\nabla K(p_1)-\nabla K(p(t))\|\le\gamma\|p_1-p(t)\|\le\gamma Mhであるから

∥q(t+h)−q1∥≤∥q(t+h)−q(t)−h∇K(p(t))∥+h∥∇K(p1)−∇K(p(t))∥≤32γMh2,∥p(t+h)−p1∥≤γM2h2\|q(t+h)-q_1\|\le\|q(t+h)-q(t)-h\nabla K(p(t))\|+h\|\nabla K(p_1)-\nabla K(p(t))\|\le\frac32\gamma Mh^2,\qquad\|p(t+h)-p_1\|\le\frac{\gamma M}2h^2

であり、∥z(t+h)−ΦhVK(z(t))∥≤∥q(t+h)−q1∥+∥p(t+h)−p1∥≤2γMh2\|z(t+h)-\Phi^{VK}_h(z(t))\|\le\|q(t+h)-q_1\|+\|p(t+h)-p_1\|\le2\gamma Mh^2である。Φh=ΦhKV\Phi_h=\Phi^{KV}_hの場合、q1:=q(t)+h∇K(p(t))q_1:=q(t)+h\nabla K(p(t))とp1:=p(t)−h∇V(q1)p_1:=p(t)-h\nabla V(q_1)は同様にq1∈Bˉ(q(t),ρ)q_1\in\bar B(q(t),\rho)、∥p1−p(t)∥≤hM≤ρ\|p_1-p(t)\|\le hM\le\rhoを満たし、z(t)∈WhKVz(t)\in W^{KV}_hである。∥∇V(q1)−∇V(q(t))∥≤γMh\|\nabla V(q_1)-\nabla V(q(t))\|\le\gamma Mhから∥q(t+h)−q1∥≤γMh2/2\|q(t+h)-q_1\|\le\gamma Mh^2/2、∥p(t+h)−p1∥≤3γMh2/2\|p(t+h)-p_1\|\le3\gamma Mh^2/2であり、同じ上界を得る。

(4)を示す。z′(τ)=(−sin⁡τ,−cos⁡τ)z'(\tau)=(-\sin\tau,-\cos\tau)は(p(τ),−q(τ))(p(\tau),-q(\tau))に等しく、w′(τ)=(cos⁡τ,−sin⁡τ)w'(\tau)=(\cos\tau,-\sin\tau)も同じ関係を満たす。§E4.5 定理 1.1をn=1n=1、r=1r=1でcos⁡\cosに適用すると

cos⁡h=1−h2∫01(1−s)cos⁡(sh) ds≥1−h22\cos h=1-h^2\int_0^1(1-s)\cos(sh)\,ds\ge1-\frac{h^2}2

である。(2)の証明の行列によりΦhVK(z(0))=ΦhVK(1,0)=(1−h2,−h)\Phi^{VK}_h(z(0))=\Phi^{VK}_h(1,0)=(1-h^2,-h)であるから、z(h)−ΦhVK(z(0))z(h)-\Phi^{VK}_h(z(0))の第一成分はcos⁡h−1+h2≥h2/2\cos h-1+h^2\ge h^2/2である。同様にΦhKV(w(0))=ΦhKV(0,1)=(h,1−h2)\Phi^{KV}_h(w(0))=\Phi^{KV}_h(0,1)=(h,1-h^2)であるから、w(h)−ΦhKV(w(0))w(h)-\Phi^{KV}_h(w(0))の第二成分はcos⁡h−1+h2≥h2/2\cos h-1+h^2\ge h^2/2である。C≥0C\ge0とh0>0h_0>0を取り、0<h<min⁡{h0,1,1/(2C)}0<h<\min\{h_0,1,1/(2C)\}（C=0C=0のときは0<h<min⁡{h0,1}0<h<\min\{h_0,1\}）とするとCh3<h2/2Ch^3<h^2/2であるから、最後の主張を得る。▨

定理 2.6.定義 1.2の設定で、Rd\R^dとR2d=Rd×Rd\R^{2d}=\R^d\times\R^dに Euclid ノルムを入れ、Sh ⁣:WhSV→R2dS_h\colon W^{\mathrm{SV}}_h\to\R^{2d}を Störmer–Verlet 法の一歩写像とする。

  1. 任意のh∈Rh\in\Rについて、ShS_hは symplectic 写像であり、Sh(WhSV)=W−hSVS_h(W^{\mathrm{SV}}_h)=W^{\mathrm{SV}}_{-h}かつ任意のz∈WhSVz\in W^{\mathrm{SV}}_hについてS−h(Sh(z))=zS_{-h}(S_h(z))=zである。(Sh)h∈R(S_h)_{h\in\R}は時間対称である。
  2. K∈C3(Up;R)K\in C^3(U_p;\R)、V∈C3(Uq;R)V\in C^3(U_q;\R)とし、t0<Tt_0<T、[t0,T][t_0,T]を含む開区間JJと、Hamilton 系を満たすC1C^1級写像z=(q,p) ⁣:J→Uz=(q,p)\colon J\to Uを取る。ρ>0\rho>0とM>0M>0が次を満たすとする。任意のτ∈[t0,T]\tau\in[t_0,T]について閉球Bˉ(q(τ),ρ)\bar B(q(\tau),\rho)はUqU_qに、Bˉ(p(τ),ρ)\bar B(p(\tau),\rho)はUpU_pに含まれ、これらの閉球のτ∈[t0,T]\tau\in[t_0,T]にわたる和集合をそれぞれNqN_q、NpN_pとすると、任意のx∈Nqx\in N_q、y∈Npy\in N_pとj∈{1,2}j\in\{1,2\}について ∥∇V(x)∥≤M,∥Dj(∇V)(x)∥op≤M,∥∇K(y)∥≤M,∥Dj(∇K)(y)∥op≤M\|\nabla V(x)\|\le M,\quad\|D^j(\nabla V)(x)\|_{\mathrm{op}}\le M,\quad\|\nabla K(y)\|\le M,\quad\|D^j(\nabla K)(y)\|_{\mathrm{op}}\le M が成り立つ。γ:=d M\gamma:=\sqrt d\,M、H0:=3ρ/(4M)H_0:=3\rho/(4M)、Λ:=2γ+54γ2H0+14γ3H02\Lambda:=2\gamma+\frac54\gamma^2H_0+\frac14\gamma^3H_0^2と置く。このとき、t∈[t0,T]t\in[t_0,T]、0<h≤H00<h\le H_0と∥u−z(t)∥≤ρ/4\|u-z(t)\|\le\rho/4を満たす任意のt,h,ut,h,uについて u∈WhSV,∥Sh(u)−Sh(z(t))∥≤(1+Λh)∥u−z(t)∥u\in W^{\mathrm{SV}}_h,\qquad\|S_h(u)-S_h(z(t))\|\le(1+\Lambda h)\|u-z(t)\| が成り立ち、t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{H0,T−t}0<h\le\min\{H_0,T-t\}を満たす任意のt,ht,hについて ∥z(t+h)−Sh(z(t))∥≤2γM2h3\|z(t+h)-S_h(z(t))\|\le2\gamma M^2h^3 が成り立つ。
  3. (2)の仮定の下でφ:=(eΛ(T−t0)−1)/Λ\varphi:=(e^{\Lambda(T-t_0)}-1)/\Lambdaと置く。[t0,T][t_0,T]の格子t0<t1<⋯<tN=Tt_0<t_1<\dots<t_N=Tの刻みhn:=tn+1−tnh_n:=t_{n+1}-t_nがH:=max⁡nhn≤H0H:=\max_nh_n\le H_0と2γM2φH2≤ρ/42\gamma M^2\varphi H^2\le\rho/4を満たすならば、z0:=z(t0)z_0:=z(t_0)とzn+1:=Shn(zn)z_{n+1}:=S_{h_n}(z_n)で定まるz0,…,zNz_0,\dots,z_Nはすべて定まり、 max⁡0≤n≤N∥zn−z(tn)∥≤2γM2φH2\max_{0\le n\le N}\|z_n-z(t_n)\|\le2\gamma M^2\varphi H^2 が成り立つ。

証明.(1)を示す。補題 1.4 (3)を列(V,K,V)(V,K,V)と(h/2,h,h/2)(h/2,h,h/2)に適用すると、S=ShS=S_h、W=WhSVW=W^{\mathrm{SV}}_hであり、逆順の列(V,K,V)(V,K,V)と(−h/2,−h,−h/2)(-h/2,-h,-h/2)から定まる写像と集合はS−=S−hS^-=S_{-h}、W−=W−hSVW^-=W^{\mathrm{SV}}_{-h}である。したがってShS_hは symplectic 写像であり、Sh(WhSV)=W−hSVS_h(W^{\mathrm{SV}}_h)=W^{\mathrm{SV}}_{-h}かつS−h∘Sh=id⁡WhSVS_{-h}\circ S_h=\id_{W^{\mathrm{SV}}_h}である。hhを−h-hに替えるとS−h(W−hSV)=WhSVS_{-h}(W^{\mathrm{SV}}_{-h})=W^{\mathrm{SV}}_hかつSh∘S−h=id⁡W−hSVS_h\circ S_{-h}=\id_{W^{\mathrm{SV}}_{-h}}であるから、S−hS_{-h}はW−hSVW^{\mathrm{SV}}_{-h}からWhSVW^{\mathrm{SV}}_hへの全単射であり、Sh∗=(S−h)−1=ShS^*_h=(S_{-h})^{-1}=S_hである。

(2)を示す。t∈[t0,T]t\in[t_0,T]を取る。Bˉ(q(t),ρ)\bar B(q(t),\rho)とBˉ(p(t),ρ)\bar B(p(t),\rho)は凸であり、それぞれNqN_q、NpN_pに含まれるから、補題 2.4をr=0r=0で適用すると、∇V\nabla VはBˉ(q(t),ρ)\bar B(q(t),\rho)上で、∇K\nabla KはBˉ(p(t),ρ)\bar B(p(t),\rho)上で、それぞれγ\gammaを係数とする Lipschitz 条件を満たす。0<h≤H00<h\le H_0とし、∥u−z(t)∥≤ρ/4\|u-z(t)\|\le\rho/4を満たすu=(qu,pu)u=(q^u,p^u)を取る。hM≤3ρ/4hM\le3\rho/4である。

pˉu:=pu−h2∇V(qu),q1u:=qu+h∇K(pˉu),p1u:=pˉu−h2∇V(q1u)\bar p^u:=p^u-\tfrac h2\nabla V(q^u),\qquad q_1^u:=q^u+h\nabla K(\bar p^u),\qquad p_1^u:=\bar p^u-\tfrac h2\nabla V(q_1^u)

と置き、u=z(t)u=z(t)のときの三点をpˉ\bar p、q1q_1、p1p_1と書く。∥qu−q(t)∥≤ρ/4\|q^u-q(t)\|\le\rho/4と∥pu−p(t)∥≤ρ/4\|p^u-p(t)\|\le\rho/4であるからqu∈Bˉ(q(t),ρ)q^u\in\bar B(q(t),\rho)、pu∈Bˉ(p(t),ρ)p^u\in\bar B(p(t),\rho)、u∈Uu\in Uかつ∥∇V(qu)∥≤M\|\nabla V(q^u)\|\le Mであり、∥pˉu−p(t)∥≤ρ/4+hM/2≤ρ\|\bar p^u-p(t)\|\le\rho/4+hM/2\le\rhoであるからpˉu∈Bˉ(p(t),ρ)\bar p^u\in\bar B(p(t),\rho)かつ∥∇K(pˉu)∥≤M\|\nabla K(\bar p^u)\|\le Mである。したがって∥q1u−q(t)∥≤ρ/4+hM≤ρ\|q_1^u-q(t)\|\le\rho/4+hM\le\rhoでありq1u∈Bˉ(q(t),ρ)q_1^u\in\bar B(q(t),\rho)、さらに∥p1u−p(t)∥≤ρ/4+hM≤ρ\|p_1^u-p(t)\|\le\rho/4+hM\le\rhoでありp1u∈Upp_1^u\in U_pである。各部分流の定義域の条件が成り立つからu∈WhSVu\in W^{\mathrm{SV}}_hかつSh(u)=(q1u,p1u)S_h(u)=(q_1^u,p_1^u)であり、u=z(t)u=z(t)についても同じである。e0:=∥u−z(t)∥e_0:=\|u-z(t)\|とすると、

(qu−q(t), pˉu−pˉ)=(u−z(t))−h2(0, ∇V(qu)−∇V(q(t)))(q^u-q(t),\ \bar p^u-\bar p)=(u-z(t))-\tfrac h2\bigl(0,\ \nabla V(q^u)-\nabla V(q(t))\bigr)

と Lipschitz 条件によりe1:=∥(qu−q(t),pˉu−pˉ)∥≤(1+γh/2)e0e_1:=\|(q^u-q(t),\bar p^u-\bar p)\|\le(1+\gamma h/2)e_0である。pˉu,pˉ∈Bˉ(p(t),ρ)\bar p^u,\bar p\in\bar B(p(t),\rho)であるから同様にe2:=∥(q1u−q1,pˉu−pˉ)∥≤(1+γh)e1e_2:=\|(q_1^u-q_1,\bar p^u-\bar p)\|\le(1+\gamma h)e_1であり、q1u,q1∈Bˉ(q(t),ρ)q_1^u,q_1\in\bar B(q(t),\rho)であるから∥Sh(u)−Sh(z(t))∥≤(1+γh/2)e2\|S_h(u)-S_h(z(t))\|\le(1+\gamma h/2)e_2である。0<h≤H00<h\le H_0について

(1+γh2)2(1+γh)=1+2γh+54γ2h2+14γ3h3≤1+Λh\Bigl(1+\frac{\gamma h}2\Bigr)^2(1+\gamma h)=1+2\gamma h+\frac54\gamma^2h^2+\frac14\gamma^3h^3\le1+\Lambda h

であるから、一つ目の主張を得る。

t∈[t0,T)t\in[t_0,T)と0<h≤min⁡{H0,T−t}0<h\le\min\{H_0,T-t\}を取る。上で示したとおりSh(z(t))=(q1,p1)S_h(z(t))=(q_1,p_1)である。∇K\nabla Kと∇V\nabla VはC2C^2級であるから、§E4.3 定理 1.1によりqqとppはJJ上でC3C^3級であり、τ∈J\tau\in Jについて

q′′(τ)=−D(∇K)(p(τ))∇V(q(τ)),q′′′(τ)=D2(∇K)(p(τ))[∇V(q(τ)),∇V(q(τ))]−D(∇K)(p(τ))D(∇V)(q(τ))∇K(p(τ)),p′′(τ)=−D(∇V)(q(τ))∇K(p(τ)),p′′′(τ)=−D2(∇V)(q(τ))[∇K(p(τ)),∇K(p(τ))]+D(∇V)(q(τ))D(∇K)(p(τ))∇V(q(τ))\begin{aligned} q''(\tau)&=-D(\nabla K)(p(\tau))\nabla V(q(\tau)),\\ q'''(\tau)&=D^2(\nabla K)(p(\tau))[\nabla V(q(\tau)),\nabla V(q(\tau))]-D(\nabla K)(p(\tau))D(\nabla V)(q(\tau))\nabla K(p(\tau)),\\ p''(\tau)&=-D(\nabla V)(q(\tau))\nabla K(p(\tau)),\\ p'''(\tau)&=-D^2(\nabla V)(q(\tau))[\nabla K(p(\tau)),\nabla K(p(\tau))]+D(\nabla V)(q(\tau))D(\nabla K)(p(\tau))\nabla V(q(\tau)) \end{aligned}

である。τ∈[t0,T]\tau\in[t_0,T]についてq(τ)∈Nqq(\tau)\in N_q、p(τ)∈Npp(\tau)\in N_pであるから∥q′′′(τ)∥≤2M3\|q'''(\tau)\|\le2M^3、∥p′′′(τ)∥≤2M3\|p'''(\tau)\|\le2M^3である。[t,t+h]⊆[t0,T][t,t+h]\subseteq[t_0,T]であるから、補題 2.4をm=1m=1、n=dn=d、r=2r=2でqqとppに適用すると

∥q(t+h)−q(t)−hq′(t)−h22q′′(t)∥≤γM23h3,∥p(t+h)−p(t)−hp′(t)−h22p′′(t)∥≤γM23h3\Bigl\|q(t+h)-q(t)-hq'(t)-\frac{h^2}2q''(t)\Bigr\|\le\frac{\gamma M^2}3h^3,\qquad\Bigl\|p(t+h)-p(t)-hp'(t)-\frac{h^2}2p''(t)\Bigr\|\le\frac{\gamma M^2}3h^3

である。∥pˉ−p(t)∥≤hM/2\|\bar p-p(t)\|\le hM/2であるから、補題 2.4をr=1r=1で∇K\nabla Kに適用すると∥∇K(pˉ)−∇K(p(t))−D(∇K)(p(t))(pˉ−p(t))∥≤γM2h2/8\|\nabla K(\bar p)-\nabla K(p(t))-D(\nabla K)(p(t))(\bar p-p(t))\|\le\gamma M^2h^2/8である。hD(∇K)(p(t))(pˉ−p(t))=h22q′′(t)hD(\nabla K)(p(t))(\bar p-p(t))=\frac{h^2}2q''(t)とh∇K(p(t))=hq′(t)h\nabla K(p(t))=hq'(t)により

∥q1−q(t)−hq′(t)−h22q′′(t)∥≤γM28h3\Bigl\|q_1-q(t)-hq'(t)-\frac{h^2}2q''(t)\Bigr\|\le\frac{\gamma M^2}8h^3

である。∥q1−q(t)∥≤hM\|q_1-q(t)\|\le hMであるから、補題 2.4をr=1r=1で∇V\nabla Vに適用すると∥∇V(q1)−∇V(q(t))−D(∇V)(q(t))(q1−q(t))∥≤γM2h2/2\|\nabla V(q_1)-\nabla V(q(t))-D(\nabla V)(q(t))(q_1-q(t))\|\le\gamma M^2h^2/2である。Lipschitz 条件により∥q1−q(t)−h∇K(p(t))∥=h∥∇K(pˉ)−∇K(p(t))∥≤γMh2/2\|q_1-q(t)-h\nabla K(p(t))\|=h\|\nabla K(\bar p)-\nabla K(p(t))\|\le\gamma Mh^2/2であり、∥D(∇V)(q(t))∥op≤M\|D(\nabla V)(q(t))\|_{\mathrm{op}}\le Mであるから∥D(∇V)(q(t))(q1−q(t))−hD(∇V)(q(t))∇K(p(t))∥≤γM2h2/2\|D(\nabla V)(q(t))(q_1-q(t))-hD(\nabla V)(q(t))\nabla K(p(t))\|\le\gamma M^2h^2/2である。したがってη:=∇V(q1)−∇V(q(t))−hD(∇V)(q(t))∇K(p(t))\eta:=\nabla V(q_1)-\nabla V(q(t))-hD(\nabla V)(q(t))\nabla K(p(t))は∥η∥≤γM2h2\|\eta\|\le\gamma M^2h^2を満たし、

p1=p(t)−h∇V(q(t))−h22D(∇V)(q(t))∇K(p(t))−h2η=p(t)+hp′(t)+h22p′′(t)−h2ηp_1=p(t)-h\nabla V(q(t))-\frac{h^2}2D(\nabla V)(q(t))\nabla K(p(t))-\frac h2\eta=p(t)+hp'(t)+\frac{h^2}2p''(t)-\frac h2\eta

である。以上を合わせると

∥z(t+h)−Sh(z(t))∥≤∥q(t+h)−q1∥+∥p(t+h)−p1∥≤(13+18)γM2h3+(13+12)γM2h3=3124γM2h3\|z(t+h)-S_h(z(t))\|\le\|q(t+h)-q_1\|+\|p(t+h)-p_1\|\le\Bigl(\frac13+\frac18\Bigr)\gamma M^2h^3+\Bigl(\frac13+\frac12\Bigr)\gamma M^2h^3=\frac{31}{24}\gamma M^2h^3

であり、二つ目の主張を得る。

(3)を示す。(2)の所属u∈WhSVu\in W^{\mathrm{SV}}_hにより(t,u,h)(t,u,h)は Störmer–Verlet 法の増分関数の定義域に属するから、同命題の二つの評価により、§E20.28 系 1.7の仮定はρ\rhoをρ/4\rho/4、p=2p=2、C=2γM2C=2\gamma M^2として成り立つ。γ>0\gamma>0であるからΛ>0\Lambda>0であり、φ\varphiは同系の定数である。y0=z(t0)y_0=z(t_0)、rn=0r_n=0とするとE0=0E_0=0、ε=0\varepsilon=0であり、同系の条件の左辺は2γM2φH2≤ρ/42\gamma M^2\varphi H^2\le\rho/4となるから、z0,…,zNz_0,\dots,z_Nはすべて定まり、max⁡n∥zn−z(tn)∥≤2γM2φH2\max_n\|z_n-z(t_n)\|\le2\gamma M^2\varphi H^2である。▨

3 調和振動子

命題 3.1.d=1d=1、Uq=Up=RU_q=U_p=\R、K(p)=p2/2K(p)=p^2/2、V(q)=q2/2V(q)=q^2/2とし、H(q,p)=(q2+p2)/2H(q,p)=(q^2+p^2)/2とする。R2\R^2の元を列ベクトルとみなし、Euclid ノルムを入れる。h,s,ϕ∈Rh,s,\phi\in\Rに対して

Ah:=(1−h2h−h1),Bh:=(1−h22h−h(1−h24)1−h22),Es:=(10−s1),R(ϕ):=(cos⁡ϕsin⁡ϕ−sin⁡ϕcos⁡ϕ)A_h:=\begin{pmatrix}1-h^2&h\\-h&1\end{pmatrix},\quad B_h:=\begin{pmatrix}1-\frac{h^2}2&h\\-h\bigl(1-\frac{h^2}4\bigr)&1-\frac{h^2}2\end{pmatrix},\quad E_s:=\begin{pmatrix}1&0\\-s&1\end{pmatrix},\quad R(\phi):=\begin{pmatrix}\cos\phi&\sin\phi\\-\sin\phi&\cos\phi\end{pmatrix}

と置き、Qh(q,p):=q2+p2−hqpQ_h(q,p):=q^2+p^2-hqp、QhSV(q,p):=(1−h24)q2+p2Q^{\mathrm{SV}}_h(q,p):=\bigl(1-\frac{h^2}4\bigr)q^2+p^2と置く。

  1. 任意のh∈Rh\in\Rとz∈R2z\in\R^2について、ΦhVK(z)=Ahz\Phi^{VK}_h(z)=A_hz、Sh(z)=BhzS_h(z)=B_hz、Bh=Eh/2AhEh/2−1B_h=E_{h/2}A_hE_{h/2}^{-1}であり、 Qh(Ahz)=Qh(z),QhSV(Bhz)=QhSV(z)Q_h(A_hz)=Q_h(z),\qquad Q^{\mathrm{SV}}_h(B_hz)=Q^{\mathrm{SV}}_h(z) が成り立つ。任意のz0∈R2z_0\in\R^2について、τ↦R(τ)z0\tau\mapsto R(\tau)z_0は Hamilton 系を満たす。
  2. h≠0h\ne0ならば、e2:=(0,1)e_2:=(0,1)についてH(Ahe2)≠H(e2)H(A_he_2)\ne H(e_2)かつH(Bhe2)≠H(e2)H(B_he_2)\ne H(e_2)である。
  3. 任意のh∈Rh\in\Rとz∈R2z\in\R^2について (1−∣h∣2)∥z∥2≤Qh(z)≤(1+∣h∣2)∥z∥2,(1−h24)∥z∥2≤QhSV(z)≤∥z∥2\Bigl(1-\frac{|h|}2\Bigr)\|z\|^2\le Q_h(z)\le\Bigl(1+\frac{|h|}2\Bigr)\|z\|^2,\qquad\Bigl(1-\frac{h^2}4\Bigr)\|z\|^2\le Q^{\mathrm{SV}}_h(z)\le\|z\|^2 が成り立つ。QhQ_hとQhSVQ^{\mathrm{SV}}_hは、∣h∣<2|h|<2のとき正定値、∣h∣=2|h|=2のとき半正定値であって正定値でなく、∣h∣>2|h|>2のとき不定符号である。
  4. ∣h∣<2|h|<2とし、θ(h):=2arcsin⁡(h/2)\theta(h):=2\arcsin(h/2)、ch:=1−h2/4c_h:=\sqrt{1-h^2/4}、 Lh:=(1−h/20ch),LhSV:=LhEh/2−1L_h:=\begin{pmatrix}1&-h/2\\0&c_h\end{pmatrix},\qquad L^{\mathrm{SV}}_h:=L_hE_{h/2}^{-1} と置く。任意のz∈R2z\in\R^2についてQh(z)=∥Lhz∥2Q_h(z)=\|L_hz\|^2、QhSV(z)=∥LhSVz∥2Q^{\mathrm{SV}}_h(z)=\|L^{\mathrm{SV}}_hz\|^2であり、 LhAhLh−1=LhSVBh(LhSV)−1=R(θ(h))L_hA_hL_h^{-1}=L^{\mathrm{SV}}_hB_h(L^{\mathrm{SV}}_h)^{-1}=R(\theta(h)) が成り立つ。さらに0<∣h∣≤10<|h|\le1について0≤θ(h)/h−1−h2/24≤h4/1000\le\theta(h)/h-1-h^2/24\le h^4/100が成り立つ。
  5. ∣h∣≥2|h|\ge2ならば、z,z′∈R2z,z'\in\R^2が存在して、(∥Ahnz∥)n∈N≥0(\|A_h^nz\|)_{n\in\N}と(∥Bhnz′∥)n∈N≥0(\|B_h^nz'\|)_{n\in\N}は有界でない。

証明.(1)を示す。∇V(q)=q\nabla V(q)=q、∇K(p)=p\nabla K(p)=pであるから、すべての部分流の定義域はR2\R^2であり、φsV(z)=Esz\varphi^V_s(z)=E_sz、φhK(z)=Fhz\varphi^K_h(z)=F_hz（Fh:=(1h01)F_h:=\begin{pmatrix}1&h\\0&1\end{pmatrix}）である。Ah=FhEhA_h=F_hE_hであり、EaEb=Ea+bE_aE_b=E_{a+b}により

Sh(z)=Eh/2FhEh/2z=Eh/2FhEhE−h/2z=Eh/2AhEh/2−1zS_h(z)=E_{h/2}F_hE_{h/2}z=E_{h/2}F_hE_hE_{-h/2}z=E_{h/2}A_hE_{h/2}^{-1}z

である。右辺の積を計算するとBhB_hを得る。(q1,p1):=Ah(q,p)(q_1,p_1):=A_h(q,p)とするとq1−hp1=qq_1-hp_1=qであるから

Qh(q1,p1)=q1(q1−hp1)+p12=q((1−h2)q+hp)+(p−hq)2=q2+p2−hqpQ_h(q_1,p_1)=q_1(q_1-hp_1)+p_1^2=q\bigl((1-h^2)q+hp\bigr)+(p-hq)^2=q^2+p^2-hqp

である。E−h/2(q,p)=(q,p+hq/2)E_{-h/2}(q,p)=(q,p+hq/2)であるからQh(E−h/2z)=(1−h2/4)q2+p2=QhSV(z)Q_h(E_{-h/2}z)=(1-h^2/4)q^2+p^2=Q^{\mathrm{SV}}_h(z)であり、

QhSV(Bhz)=Qh(E−h/2Eh/2AhE−h/2z)=Qh(AhE−h/2z)=Qh(E−h/2z)=QhSV(z)Q^{\mathrm{SV}}_h(B_hz)=Q_h(E_{-h/2}E_{h/2}A_hE_{-h/2}z)=Q_h(A_hE_{-h/2}z)=Q_h(E_{-h/2}z)=Q^{\mathrm{SV}}_h(z)

である。R(τ)z0=(q(τ),p(τ))R(\tau)z_0=(q(\tau),p(\tau))と置くと、ddτR(τ)=(01−10)R(τ)\frac{d}{d\tau}R(\tau)=\begin{pmatrix}0&1\\-1&0\end{pmatrix}R(\tau)によりq′=pq'=p、p′=−qp'=-qである。

(2)を示す。Ahe2=(h,1)A_he_2=(h,1)とBhe2=(h,1−h2/2)B_he_2=(h,1-h^2/2)から2H(Ahe2)=1+h22H(A_he_2)=1+h^2、2H(Bhe2)=h2+(1−h2/2)2=1+h4/42H(B_he_2)=h^2+(1-h^2/2)^2=1+h^4/4であり、h≠0h\ne0ならばいずれも2H(e2)=12H(e_2)=1に等しくない。

(3)を示す。z=(q,p)z=(q,p)について∣qp∣≤∥z∥2/2|qp|\le\|z\|^2/2であるからQhQ_hの評価を得る。QhSV(z)−(1−h2/4)∥z∥2=h24p2≥0Q^{\mathrm{SV}}_h(z)-(1-h^2/4)\|z\|^2=\frac{h^2}4p^2\ge0と∥z∥2−QhSV(z)=h24q2≥0\|z\|^2-Q^{\mathrm{SV}}_h(z)=\frac{h^2}4q^2\ge0からQhSVQ^{\mathrm{SV}}_hの評価を得る。∣h∣<2|h|<2ならば両者の下界の係数は正であるから、二つの二次形式は正定値である。h=±2h=\pm2のときQh(q,p)=(q∓p)2Q_h(q,p)=(q\mp p)^2、QhSV(q,p)=p2Q^{\mathrm{SV}}_h(q,p)=p^2は非負であり、それぞれ(1,±1)(1,\pm1)、(1,0)(1,0)で00となる。∣h∣>2|h|>2のときQh(1,1)=2−hQ_h(1,1)=2-hとQh(1,−1)=2+hQ_h(1,-1)=2+hは異符号であり、QhSV(1,0)=1−h2/4<0<1=QhSV(0,1)Q^{\mathrm{SV}}_h(1,0)=1-h^2/4<0<1=Q^{\mathrm{SV}}_h(0,1)である。

(4)を示す。∥Lhz∥2=(q−hp/2)2+(1−h2/4)p2=Qh(z)\|L_hz\|^2=(q-hp/2)^2+(1-h^2/4)p^2=Q_h(z)であり、(1)の証明により∥LhSVz∥2=Qh(E−h/2z)=QhSV(z)\|L^{\mathrm{SV}}_hz\|^2=Q_h(E_{-h/2}z)=Q^{\mathrm{SV}}_h(z)である。

LhAhLh−1=(1−h22h2−chhch)(1h2ch01ch)=(1−h22chh−chh1−h22)L_hA_hL_h^{-1}=\begin{pmatrix}1-\frac{h^2}2&\frac h2\\-c_hh&c_h\end{pmatrix}\begin{pmatrix}1&\frac h{2c_h}\\0&\frac1{c_h}\end{pmatrix}=\begin{pmatrix}1-\frac{h^2}2&c_hh\\-c_hh&1-\frac{h^2}2\end{pmatrix}

である。ここで(1,2)(1,2)成分はh2ch(2−h22)=hchch2\frac h{2c_h}\bigl(2-\frac{h^2}2\bigr)=\frac{h}{c_h}c_h^2を用いた。∣θ(h)/2∣<π/2|\theta(h)/2|<\pi/2であるからcos⁡(θ(h)/2)=ch\cos(\theta(h)/2)=c_hであり、cos⁡θ(h)=1−2sin⁡2(θ(h)/2)=1−h2/2\cos\theta(h)=1-2\sin^2(\theta(h)/2)=1-h^2/2、sin⁡θ(h)=2sin⁡(θ(h)/2)cos⁡(θ(h)/2)=hch\sin\theta(h)=2\sin(\theta(h)/2)\cos(\theta(h)/2)=hc_hである。したがってLhAhLh−1=R(θ(h))L_hA_hL_h^{-1}=R(\theta(h))であり、

LhSVBh(LhSV)−1=LhEh/2−1(Eh/2AhEh/2−1)Eh/2Lh−1=LhAhLh−1L^{\mathrm{SV}}_hB_h(L^{\mathrm{SV}}_h)^{-1}=L_hE_{h/2}^{-1}\bigl(E_{h/2}A_hE_{h/2}^{-1}\bigr)E_{h/2}L_h^{-1}=L_hA_hL_h^{-1}

である。θ(h)/h\theta(h)/hの評価の証明は演習とする（問題 5.1）。

(5)を示す。∣h∣=2|h|=2のとき、h2=4h^2=4であるから

N:=Ah+I2=(−2h−h2),N2=(4−h2−2h+2h2h−2h4−h2)=0N:=A_h+I_2=\begin{pmatrix}-2&h\\-h&2\end{pmatrix},\qquad N^2=\begin{pmatrix}4-h^2&-2h+2h\\2h-2h&4-h^2\end{pmatrix}=0

であり、(1,2)(1,2)成分hhが00でないからN≠0N\ne0である。Ahn=(−1)n(I2−N)n=(−1)n(I2−nN)A_h^n=(-1)^n(I_2-N)^n=(-1)^n(I_2-nN)であるから、Nz≠0Nz\ne0を満たすzzについて∥Ahnz∥≥n∥Nz∥−∥z∥\|A_h^nz\|\ge n\|Nz\|-\|z\|である。∣h∣>2|h|>2のとき、tr⁡Ah=2−h2\operatorname{tr}A_h=2-h^2、det⁡Ah=1\det A_h=1であるから、AhA_hの固有多項式はλ2−(2−h2)λ+1\lambda^2-(2-h^2)\lambda+1である。その判別式(2−h2)2−4(2-h^2)^2-4は正であり、二つの実根の積は11、和は2−h2<−22-h^2<-2である。二つの根の絶対値がともに11以下ならば、積が11であるから両者は±1\pm1であり、和は−2-2以上となるので、絶対値が11より大きい根λ\lambdaが存在する。その固有ベクトルzzについて∥Ahnz∥=∣λ∣n∥z∥\|A_h^nz\|=|\lambda|^n\|z\|である。いずれの場合もz′:=Eh/2zz':=E_{h/2}zと置くとBhnz′=Eh/2AhnzB_h^nz'=E_{h/2}A_h^nzであり、∥Ahnz∥≤∥Eh/2−1∥op∥Bhnz′∥\|A_h^nz\|\le\|E_{h/2}^{-1}\|_{\mathrm{op}}\|B_h^nz'\|であるから、(∥Bhnz′∥)n(\|B_h^nz'\|)_nも有界でない。▨

例 3.2.命題 3.1の設定でh=1/10h=1/10、z0=(1,0)z_0=(1,0)とし、前進 Euler 法、ΦhVK\Phi^{VK}_h、ShS_hのそれぞれでz0z_0からnn歩進めた値をznz_nと書く。厳密解R(τ)z0R(\tau)z_0は∥R(τ)z0∥=1\|R(\tau)z_0\|=1とH(R(τ)z0)=1/2H(R(\tau)z_0)=1/2を満たし、時刻がhh進むごとに角hhだけ回転する。

前進 Euler 法の一歩はz↦z+h(p,−q)z\mapsto z+h(p,-q)であり、その行列は

(1h−h1)=1+h2 R(arctan⁡h)\begin{pmatrix}1&h\\-h&1\end{pmatrix}=\sqrt{1+h^2}\,R(\arctan h)

である。したがって∥zn∥=(1+h2)n/2\|z_n\|=(1+h^2)^{n/2}、H(zn)=(1+h2)n/2H(z_n)=(1+h^2)^n/2であり、znz_nの角は一歩ごとにarctan⁡h=0.0996687…\arctan h=0.0996687\ldotsだけ進む。n=1000n=1000（時刻100100）では∥zn∥=1.01500=144.77…\|z_n\|=1.01^{500}=144.77\ldotsであり、角のずれはn(arctan⁡h−h)=−0.3313…n(\arctan h-h)=-0.3313\ldotsである。

symplectic Euler 法ΦhVK\Phi^{VK}_hについては、命題 3.1 (4)によりLhzn=R(nθ(h))Lhz0L_hz_n=R(n\theta(h))L_hz_0であり、∥Lhzn∥2=Qh(z0)=1\|L_hz_n\|^2=Q_h(z_0)=1はすべてのnnで変わらない。命題 3.1 (3)により、すべてのn∈N≥0n\in\Nについて

11.05≤∥zn∥2≤10.95,0.47619…≤H(zn)≤0.52631…\frac1{1.05}\le\|z_n\|^2\le\frac1{0.95},\qquad0.47619\ldots\le H(z_n)\le0.52631\ldots

である。(Lhzn)n(L_hz_n)_nは一歩ごとに角θ(1/10)=0.1000417…\theta(1/10)=0.1000417\ldotsだけ回転し、n=1000n=1000での角のずれはn(θ(h)−h)=0.04171…n(\theta(h)-h)=0.04171\ldotsである。この値は命題 3.1 (4)の評価から従う区間[nh3/24, n(h3/24+h5/100)]=[0.041666…,0.041766…][nh^3/24,\ n(h^3/24+h^5/100)]=[0.041666\ldots,0.041766\ldots]に属する。

Störmer–Verlet 法ShS_hについては、命題 3.1 (1)によりQhSV(zn)=QhSV(z0)=1−h2/4=0.9975Q^{\mathrm{SV}}_h(z_n)=Q^{\mathrm{SV}}_h(z_0)=1-h^2/4=0.9975であり、命題 3.1 (3)によりすべてのnnについて0.9975≤∥zn∥2≤10.9975\le\|z_n\|^2\le1、すなわち∣H(zn)−H(z0)∣≤0.00125|H(z_n)-H(z_0)|\le0.00125である。(LhSVzn)n(L^{\mathrm{SV}}_hz_n)_nはΦhVK\Phi^{VK}_hの場合と同じ角θ(h)\theta(h)で回転する。

一般の0<∣h∣<20<|h|<2とz0∈R2z_0\in\R^2についても、命題 3.1 (1)と命題 3.1 (3)により、すべてのn∈N≥0n\in\Nで

Qh(z0)1+∣h∣/2≤2H((ΦhVK)n(z0))≤Qh(z0)1−∣h∣/2,QhSV(z0)≤2H(Shn(z0))≤QhSV(z0)1−h2/4\frac{Q_h(z_0)}{1+|h|/2}\le2H\bigl((\Phi^{VK}_h)^n(z_0)\bigr)\le\frac{Q_h(z_0)}{1-|h|/2},\qquad Q^{\mathrm{SV}}_h(z_0)\le2H(S_h^n(z_0))\le\frac{Q^{\mathrm{SV}}_h(z_0)}{1-h^2/4}

が成り立つ。HHの値はnnによらない幅に収まるが、命題 3.1 (2)によりHHはAhA_hでもBhB_hでも保存されない。

4 振り子と状態に依存する刻み

例 4.1.d=1d=1、Uq=Up=RU_q=U_p=\R、K(p)=p2/2K(p)=p^2/2、V(q)=1−cos⁡qV(q)=1-\cos qとし、H(q,p)=p2/2+1−cos⁡qH(q,p)=p^2/2+1-\cos qとする。すべての部分流の定義域はR2\R^2であり、Störmer–Verlet 法の一歩は

pˉ=p−h2sin⁡q,q1=q+hpˉ,p1=pˉ−h2sin⁡q1\bar p=p-\tfrac h2\sin q,\qquad q_1=q+h\bar p,\qquad p_1=\bar p-\tfrac h2\sin q_1

である。初期値z0=(2,0)z_0=(2,0)（H(z0)=1−cos⁡2=1.4161468…H(z_0)=1-\cos2=1.4161468\ldots）から、次の二通りに近似値を計算した。

  • 固定刻み:zn+1=S1/10(zn)z_{n+1}=S_{1/10}(z_n)、tn=n/10t_n=n/10。
  • 状態に依存する刻み:η(p):=1/51+p2\eta(p):=\frac{1/5}{1+p^2}と置き、zn=(qn,pn)z_n=(q_n,p_n)に対してzn+1=Sη(pn)(zn)z_{n+1}=S_{\eta(p_n)}(z_n)、tn+1=tn+η(pn)t_{n+1}=t_n+\eta(p_n)、t0=0t_0=0。

補題 1.3により厳密解に沿ってHHはH(z0)H(z_0)に等しい。IEEE 754 倍精度の浮動小数点演算で計算し、tn≥Tt_n\ge Tとなる最初のnnまでのmax⁡m∣H(zm)−H(z0)∣\max_m|H(z_m)-H(z_0)|を有効数字4桁に丸めると次のとおりである。

TT 固定刻みの歩数 固定刻みの最大値 状態依存刻みの歩数 tnt_n 状態依存刻みの最大値 状態依存刻みのH(zn)−H(z0)H(z_n)-H(z_0)
1010 100100 2.705×10−32.705\times10^{-3} 107107 10.0410.04 1.129×10−31.129\times10^{-3} −9.487×10−4-9.487\times10^{-4}
100100 10001000 2.706×10−32.706\times10^{-3} 11021102 100.18100.18 3.219×10−33.219\times10^{-3} −2.273×10−3-2.273\times10^{-3}
10001000 1000010000 2.706×10−32.706\times10^{-3} 1097310973 1000.011000.01 2.518×10−22.518\times10^{-2} −2.505×10−2-2.505\times10^{-2}
1000010000 100000100000 2.706×10−32.706\times10^{-3} 103872103872 10000.0410000.04 4.782×10−14.782\times10^{-1} −4.782×10−1-4.782\times10^{-1}

この入力と演算では、固定刻みの∣H(zn)−H(z0)∣|H(z_n)-H(z_0)|は時刻10410^4まで2.706×10−32.706\times10^{-3}以下にとどまった。状態依存刻みでは、表の四つの時刻でH(zn)−H(z0)H(z_n)-H(z_0)は負であり、その絶対値は時刻とともに増大した。表の四つの時刻での平均の刻みtn/nt_n/nは0.0900.090から0.0970.097の間にあり、固定刻み1/101/10より小さい。

状態依存刻みの一歩は写像F(z):=Sη(p)(z)F(z):=S_{\eta(p)}(z)である。G(h,z):=Sh(z)G(h,z):=S_h(z)はR×R2\R\times\R^2上でC1C^1級であるから、§E4.3 定理 1.1により

DF(z)=DzG(η(p),z)+∂hG(η(p),z) cT,c:=(0,η′(p))TDF(z)=D_zG(\eta(p),z)+\partial_hG(\eta(p),z)\,c^{\mathsf T},\qquad c:=(0,\eta'(p))^{\mathsf T}

である。第一項をA=(aij)A=(a_{ij})、∂hG(η(p),z)\partial_hG(\eta(p),z)をb=(bi)b=(b_i)と書く。2×22\times2行列XXについてXTJX=(det⁡X)JX^{\mathsf T}JX=(\det X)Jであるから、開集合W⊆R2W\subseteq\R^2上でFFが symplectic であることは、WWの各点でdet⁡DF=1\det DF=1であることと同値である。定理 2.6 (1)によりATJA=JA^{\mathsf T}JA=Jであるからdet⁡A=1\det A=1であり、2×22\times2行列の行列式を成分で展開すると

det⁡DF(z)=det⁡A+cT(a22−a12−a21a11)b=1+η′(p)(a11b2−a21b1)\det DF(z)=\det A+c^{\mathsf T}\begin{pmatrix}a_{22}&-a_{12}\\-a_{21}&a_{11}\end{pmatrix}b=1+\eta'(p)\bigl(a_{11}b_2-a_{21}b_1\bigr)

である。G(h,z)=(q1,p1)G(h,z)=(q_1,p_1)の成分の偏導関数は

∂q1∂q=1−h22cos⁡q,∂p1∂q=−h2cos⁡q−h2cos⁡q1∂q1∂q,∂q1∂h=p−hsin⁡q,∂p1∂h=−12sin⁡q−12sin⁡q1−h2cos⁡q1∂q1∂h\frac{\partial q_1}{\partial q}=1-\frac{h^2}2\cos q,\qquad\frac{\partial p_1}{\partial q}=-\frac h2\cos q-\frac h2\cos q_1\frac{\partial q_1}{\partial q},\qquad\frac{\partial q_1}{\partial h}=p-h\sin q,\qquad\frac{\partial p_1}{\partial h}=-\frac12\sin q-\frac12\sin q_1-\frac h2\cos q_1\frac{\partial q_1}{\partial h}

である。z=(π/2,1)z=(\pi/2,1)ではη(1)=1/10\eta(1)=1/10、η′(1)=−1/10\eta'(1)=-1/10、q1=π/2+19/200q_1=\pi/2+19/200であり、

a11=1,a21=120sin⁡19200,b1=910,b2=−12−12cos⁡19200+9200sin⁡19200a_{11}=1,\quad a_{21}=\frac1{20}\sin\frac{19}{200},\quad b_1=\frac9{10},\quad b_2=-\frac12-\frac12\cos\frac{19}{200}+\frac9{200}\sin\frac{19}{200}

であるから

det⁡DF(π2,1)=1+120(1+cos⁡19200)=1.0997745…\det DF\Bigl(\frac\pi2,1\Bigr)=1+\frac1{20}\Bigl(1+\cos\frac{19}{200}\Bigr)=1.0997745\ldots

である。DFDFは連続であるからdet⁡DF\det DFは(π/2,1)(\pi/2,1)の近傍で11と異なり、FFはその近傍で symplectic 写像でない。

5 演習

問題 5.1.命題 3.1 (4)の評価0≤θ(h)/h−1−h2/24≤h4/1000\le\theta(h)/h-1-h^2/24\le h^4/100（0<∣h∣≤10<|h|\le1）の証明を完成させよ。

解答.

θ\thetaは奇関数であるからθ(h)/h\theta(h)/hは偶関数であり、0<h≤10<h\le1としてよい。x:=h/2∈(0,1/2]x:=h/2\in(0,1/2]と置くとθ(h)/h=arcsin⁡x/x\theta(h)/h=\arcsin x/x、h2/24=x2/6h^2/24=x^2/6である。ψ(u):=(1−u)−1/2\psi(u):=(1-u)^{-1/2}は(−∞,1)(-\infty,1)上でC2C^2級であり、ψ′′(v)=34(1−v)−5/2\psi''(v)=\frac34(1-v)^{-5/2}である。§E4.5 定理 1.1をn=1n=1、r=1r=1、a=0a=0でψ\psiに適用すると

ψ(u)=1+u2+u2∫01(1−τ)ψ′′(τu) dτ\psi(u)=1+\frac u2+u^2\int_0^1(1-\tau)\psi''(\tau u)\,d\tau

である。0≤u≤1/40\le u\le1/4ならば0≤τu≤1/40\le\tau u\le1/4であり0≤ψ′′(τu)≤34(43)5/2=(43)3/2<1.540\le\psi''(\tau u)\le\frac34\bigl(\frac43\bigr)^{5/2}=\bigl(\frac43\bigr)^{3/2}<1.54であるから、0≤ψ(u)−1−u/2≤0.77u20\le\psi(u)-1-u/2\le0.77u^2である。0≤s≤x0\le s\le xについてu=s2≤1/4u=s^2\le1/4として

0≤11−s2−1−s22≤0.77s40\le\frac1{\sqrt{1-s^2}}-1-\frac{s^2}2\le0.77s^4

であり、arcsin⁡x=∫0x(1−s2)−1/2 ds\arcsin x=\int_0^x(1-s^2)^{-1/2}\,dsを[0,x][0,x]で積分すると0≤arcsin⁡x−x−x3/6≤0.154x50\le\arcsin x-x-x^3/6\le0.154x^5である。両辺をxxで割ると

0≤θ(h)h−1−h224=arcsin⁡x−x−x3/6x≤0.154x4=0.15416h4<h41000\le\frac{\theta(h)}h-1-\frac{h^2}{24}=\frac{\arcsin x-x-x^3/6}x\le0.154x^4=\frac{0.154}{16}h^4<\frac{h^4}{100}

である。▨

問題 5.2.命題 2.2の証明を完成させよ。

解答.

命題 2.2 (1)を示す。Up=RdU_p=\R^dであるから、任意の(q,p)∈U(q,p)\in Uについてp−a∇V(q)∈Upp-a\nabla V(q)\in U_pであり、WaV=UW^V_a=Uである。φaV(φbV(q,p))=(q,p−b∇V(q)−a∇V(q))=φa+bV(q,p)\varphi^V_a(\varphi^V_b(q,p))=(q,p-b\nabla V(q)-a\nabla V(q))=\varphi^V_{a+b}(q,p)である。

命題 2.2 (2)を示す。z′∈Uz'\in Uとw′:=φ−h/2V(z′)w':=\varphi^V_{-h/2}(z')について、z′∈WhSVz'\in W^{\mathrm{SV}}_hであることはφh/2V(z′)∈WhK\varphi^V_{h/2}(z')\in W^K_hと同値である。実際、キックの定義域はUUであり、ドリフトの像はUUに属するから、定義域の条件はドリフトの条件だけである。命題 2.2 (1)によりφhV(w′)=φh/2V(z′)\varphi^V_h(w')=\varphi^V_{h/2}(z')であるから、w′∈WhVKw'\in W^{VK}_hであることもφh/2V(z′)∈WhK\varphi^V_{h/2}(z')\in W^K_hと同値である。この条件の下で、命題 2.2 (1)により

Sh(z′)=φh/2V(φhK(φhV(w′)))=φh/2V(ΦhVK(w′)),φ−h/2V(Sh(z′))=ΦhVK(w′)S_h(z')=\varphi^V_{h/2}\bigl(\varphi^K_h(\varphi^V_h(w'))\bigr)=\varphi^V_{h/2}\bigl(\Phi^{VK}_h(w')\bigr),\qquad\varphi^V_{-h/2}(S_h(z'))=\Phi^{VK}_h(w')

である。0≤k≤n0\le k\le nについて、定まる限りzk:=Shk(z)z_k:=S_h^k(z)、wk:=(ΦhVK)k(w)w_k:=(\Phi^{VK}_h)^k(w)と置く。zkz_kとwkw_kが定まりwk=φ−h/2V(zk)w_k=\varphi^V_{-h/2}(z_k)であるとすると、zk∈WhSVz_k\in W^{\mathrm{SV}}_hとwk∈WhVKw_k\in W^{VK}_hは同値であり、そのときwk+1=φ−h/2V(zk+1)w_{k+1}=\varphi^V_{-h/2}(z_{k+1})である。k=0k=0でw0=φ−h/2V(z0)w_0=\varphi^V_{-h/2}(z_0)であるから、kkについての帰納法により同値性を得る。そのときzn=φh/2V(wn)z_n=\varphi^V_{h/2}(w_n)である。▨

前提記事