数理解析研究所講究録 1271 巻 2002 年 240-247240
非線形波動系に対するシンプレクティック数値解法とその拡張
日本原子力研究所

佐々成正(Narimasa Sasa)

1 はじめに常微分方程式や偏微分方程式の初期値問題を数値的に解く必要のある場合, 広く、用いられてぃる手法として, (4 次の古典)Runge-Kutta 法等の汎用パッケージを使った数値シミュレーションが挙げられる. あまり精密な結果を要求せず, 手っ取り早く何らかの結果を出すという目的においては, このような手法はある意味において優れていると言える. ただし, このような汎用的な数値解法は元の微分方程式の性質を何も反映してぃ

ないことから, 本来は現れないはずの計算結果が出ることがしばしばある. 以下では, その具体例を見てぃくことにしよう [1].
いま, Hamilton 系 (エネルギー保存系) の簡単な例として調和振動子,$\dot{q}=p$,$\dot{p}=-q$

(1)
,

(2)

を考える. Hamilton 系であるから, 当然この系ではエネルギー $H=(q^{2}+p^{2})/2$ が保存してぃなければならない.

まず, この系に前進差分 (Euler 法) を適用してみると,$q’=q+(\Delta t)p$ ,(3)$p’=p-(\Delta t)q$ ,(4)というスキームが得られる (’ は次の時間ステップにおける値を意味する).

$\cdot$

このときエネルギーは(5)$H’=(1\dotplus\Delta t^{2})H$となる. すなわち, 1 ステップ進むごとにエネルギーが $(1+\Delta t^{2})$ 倍されてぃく.

$\Delta t$

をどんなに小さく取っても必ずエネルギーは増大していくことから, 本来エネルギーが保存すべき系に対する数値計算スキームとしては適当でないことがわかる. 次に, (4 次の古典)Runge-Kutta 法を適用してみよう. 方程式 (1), (2) は線形方程式であることから時間発展の式を明示的に書くことができて,(6)$q’=(1-\Delta t^{2}/2+\Delta t^{4}/24)q+(\Delta t-\Delta t^{3}/6)p$$p’=-(\Delta t-\Delta t^{3}/6)q+(1-\Delta t^{2}/2+\Delta t^{4}/24)p$,(7)というスキームになる. このときエネルギーは(8)$H’=(1-\Delta t^{6}/72+\Delta t^{8}/576)H$で与えられ, 今度は 1 ステップ進むごとにエネルギー減少していくことになる. このように、 本来エネルギーが保存すべき系に対して, たとえが小さくても計算が長時間に及べば誤 は大きくなるので, (4 次の古典)Runge-Kutta 法も Hamilton 系に対する長時間の数値計算手法としては適当で無いことが分かる.
一方, 元の微分方程式が持っている性質を十分考慮してそれを数値計算スキームに反映させ得る場合がある. もちろん, 全ての性質をそのまま数値計算スキームに載せることはできないので, いくっかの性質を選ん$\dot{\text{差}}$$\Delta t$

でそれを反映させることになる. 例えば, 方程式 (1), (2) に対してエネルギー保存に着目したとする. いま,

陰的中点法を適用すると$q’=q+\Delta t(p’+p)/2$ ,(9)$p’=p-\Delta t(q’+q)/2$ ,

(10)

241
$\mathrm{g}r_{X}$$\hslash\grave{\grave{>}},$$\sim-\emptyset\ovalbox{\tt\small REJECT}_{\square }^{\mathrm{A}}[] \mathrm{I}$(11)$H’=H$となり, エネルギーが厳密に保存していることがわかる. また, 方程式 (1), (2) に対してハミルトン系における Symplectic 2 次形式が保存するように作った数値計算スキーム (Symplectic 数値解法) は,$q’=q+(\Delta t)p’$ ,(12)$p’=p-(\Delta t)q$ ,(13)で与えられる. このスキームは Symplectic 2 次形式が保存するだけではなく, 長時間の数値計算においてエ( I ま数値解法の次数) に比例する範囲に留まる性質があるため, エネルギーが増大や減少していくという心配は全くない.
この例からもわかるように, 方程式 (1), (2) のような非常に簡単な系に対してすら, 微分方程式が本来持っネルギーの誤差が$(\Delta t)^{n}$$n$ている性質を十分に考慮しなければ数値計算の結果が信用できるものとはなり得ないことがわかる. そこで,元の微分方程式が本来持っているいくつかの性質に着目し, それらを反映させる数値計算手法が重要となる.
これらを総称して Geometric Integrator と呼ぶことがある. 本研究では, Geometric Integrator の中の代表格である Symplectic 数値解法を非線形波動方程式系に適用する問題について考察をおこなう. また, Symplectic

数値解法はハミルトン系専用の数値解法であるが, ハミルトン系から少しずれた場合にもその計算手順が有効であることを示す.

2Symplectic 数値解法
Symplectic 数値解法は Hamilton 力学系に対して非常に有効な数値計算手法であることが知られいる [1-5].
本節では, この方法を非線形波動方程式系に適用することについてまとめてみたい [6]. ここでは例として, (複素) 従属変数 に対する次のクラスの非線形波動方程式,$u$$\mathrm{i}\frac{\partial u}{\partial t}+L[u]+W(|u|^{2})u=0$

について考察する.

(14),

ここで $L[u]$ は線形分散項で,$L[u]=\alpha_{2^{\frac{\partial^{2}u}{\partial x^{2}}}}+\cdot\alpha_{3^{\frac{\partial^{3}u}{\partial x^{3}}}}+\alpha_{4^{\frac{\partial^{4}u}{\partial x^{4}}}}+\cdots$($\alpha_{j}$

(15)

は実数)

と与えられるとする. また $W(z)$ は非線型項を表すある関数であるとする.

Symplectic 解法の一般論から, 任意の Hamilton 系に対して常に陰的な Symplectic 解法が構成できるこがとが知られてぃる $[3,4]$ . 一方, 2、ミルトユアンが $H=H_{1}+H_{2}$ の様に 2 っの和で書かれていて,共に単独では可積分 (解が簡単に書ける) であるとき, この系に対して陽的な Symplectic 解法を構成するこ$H_{1},$とが可能である.
ルトニアンとしてる.

1 次の陽的解法はまず$\Delta t$$H_{1}$

をハミルトニアンとして$\Delta t$

だけ時間発展させ, 次に$H_{2}$$H_{2}$を$’\backslash$ミ

だけ時間発展させたものの合成 (あるいはその逆) として時間 1 ステツプが定義されこれは互いに非可換な作用素 $A,$ $B$ の指数関数に関する関係式(16)$\exp[(A+B)\Delta t]=\exp[A\Delta t]\exp[B\Delta t]+O(\Delta t^{2})$がその基礎になっている. ここで $A,$ $B$ はそれぞれ$H_{1},$$A:=\{\cdot, H_{1}\},$$H_{2}$との Poisson 括弧をつくる微分作用素$B:=\{\cdot, H_{2}\}$を表す. (16) の左辺の $\exp[(A+B)\Delta t]$ は真の解の時間発展, 右辺の $\exp[A\Delta t]\exp[B\Delta t]$ はであることから,るフローの合成を表し, その差が$H_{1}$と$H_{2}$によ$O(\Delta t^{2})$$S_{1}(\Delta t):=\exp[A\Delta t]\exp[B\Delta t]$

(17)

242
が 1 次の解法であることが保証される.

またこの定義自身, 任意のハミルトニアンの解, およひその合成はsymplectic であることから, 上に述べた 1 次の解法が Symplectic 解法となることの根拠も与えている.
指数関数表記を使えば 2 次の陽解法は(18)$S_{2}(\Delta t):=\exp[A\Delta t/2]\exp[B\Delta t]\exp[A\Delta t/2]$

で与えられる.

この意味はるという意味である.

$H_{1}$

で $\Delta t/2,$

で$H_{2}$$\Delta t$

,

さらに

この

で $\Delta t/2$ 時間発展させて 1 ステツプを構成す$H_{1}$

この 2 次の陽解法は Symplectic 解法となるのみならず,

時間反転対称性を持つ解法$S_{2}(\Delta t)S_{2}(-\Delta t)=S_{2}(-\Delta t)S_{2}(\Delta t)=identity$となることも注意しておく.

4 次以上の時間反転対称な Symplectic 解法は上の 2 次の解法の対称的な合成で得られることが知られている [7]. 例えば 4 次の解法は 3 つの 2 次の解法の合成(19)$S_{4}(\Delta t):=S_{2}(x\Delta t)S_{2}((1-2x)\Delta t)S_{2}(x\Delta t)$

ここで $x=1/(2-2^{1/3})$ である.

で得られる.

ここでの結果を用いて, 非線形波動方程式 (14) に対する陽的な Symplectic 解法を構成する. まず, 方程式 (14) をに対応する部分,$H_{1}$(加)$\ovalbox{\tt\small REJECT}.\frac{\theta u}{\partial t}+\alpha_{2}\frac{\partial^{2}u}{\partial x^{2}}+\cdot\alpha_{3}\frac{\partial^{3}u}{\partial x^{3}}+\alpha_{4}\frac{\partial^{4}u}{\partial x^{4}}+\cdots=0$と$H_{2}$に対応する部分,に分ける. 方程式 (20) で.
$\Delta t$$\frac{\partial u}{\partial t}+W(|u|^{2})u=0$

(21)

,

だけ時間発展させることが $\exp[A\Delta t]$ に対応し, 方程式 (21) でさせることが $\exp[B\Delta t]$ に対応している. この対応関係と (17)-$\Delta t$

だけ時間発展(19) 等の公式を使って、 方程式 (14) に対する陽的な Symplectic 解法を構成することができる.

3Splitting スキーム
第 1 節で述べたように、 Symplectic 数値解法は Hamilton 系専用の数値計算スキームであるから, 非 Hamil-

ton 系には適用できない. しかし Hamilton 系の非線形波動方程式に簡単な散逸項を付加して、 その振舞を調べるというようなことはしばしば行われており、微小な散逸項を付け加えただけで数値計算スキームまで変更しなけれがならないとなると不便である.

そこで, 散逸を含んだ非線形波動方程式に前節で説明した陽的な Symplectic 数値解法の計算手順を適用
してみる.

このスキームは Splitting スキームと呼ばれている. 簡単な例を用いて考察してみよう. 非線形

Schr\"odinger 方程式に線形散逸項を加えた系,$\mathrm{i}\frac{\partial u}{\partial t}+\frac{\partial^{2}u}{\partial x^{2}}+2|u|^{2}u+\gamma u=0$($\gamma$: 非負の実数),(22)について考える. まず, Splitting スキームを適用するために方程式 (22) を 2 つの方程式,,.亙
笛 $+2|u|^{2}u+\gamma u=0$ ,$+ \frac{\partial^{2}u}{\partial x^{2}}=0$に分ける. 前節での議論を参考にして, 方程式 (23) で

程式 (24) で$\Delta t$$\Delta t$

(23)
(24)

だけ時間発展させることが $\exp[A\Delta t]$ に対応し, 方だけ時間発展させることが $\exp[B\Delta t]$ に対応すると考える. このとき, 例えば公式 (17) を使うと 1 次の解法が構成でき、 (18) を使うと 2 次の解法が構成できる. さらに, 前節と同じ計算手順を適用すれ

243
ば, さらに高次の計算スキームを構成することが可能である. しがし, 二般に 3 次以上のスキームの場合、 負の時間ステップが含まれる場合があって、 それが原因で数値的不安定性が起こる場合がある [8].
しかし, 方程式 (22) に対して我々が Splitting 法を適用した限りにおいてはこのような不安定性は現れな
かった.

具体的な数値計算結果を見てみょう. Fig.l は方程式 (22) に対して、 初期条件として 1 パルスの波を入れて時間発展を追ったものである. 空間領域は $-25\leq x\leq 25$ で空間メッシュを$u=\cosh(x)\exp(0.8\mathrm{i}x)$にとり, 周期境界条件とする. 空間方向の計算は FFT を用いておこなう. また, 時間 $0\leq t\leq 200$
に対して数値計算をおこない、 時間メッシュは $\Delta t=0.1$ にとる. $\gamma=0.01$ とした. 数値計算は 4 次の Splitting
法 (19) で計算をおこなった.
$\Delta x=0.195$$|\mathrm{u}|$Fig.l

1 パルスからの時間発展の様子

数値計算の結果を見ると、 線形の散逸項のためにパルスの振幅が段々と減ってくることが分がる。 実は方程式 (22) で時間発展させると$I_{1}(t)$$=$$\int|u|^{2}dx$$I_{2}(t)$

$=$

2

(25)
(26)

$\int{\rm Im}(u^{*}\frac{\partial u}{\partial x})dx$という量が, $I_{j}(t)=I_{j}(0)\exp(-\gamma t)(j=1,2)$ に従って減少することが簡単な計算から示せる。 そこで, (4 次の) Splitting 法と (4 次の) Runge-Kutta 法で時間発展させた時のの値をプロットしたのが Fig 2 で,各時間で $\tilde{I}_{10}(t)=\tilde{I}_{1}(t)-I_{1}(0)\exp(-\gamma t)$ をプロットしたものが Fig 3 である. (ここで、は (25) の数値計算における値とする). (4 次の) Runge-Kutta 法では $\Delta t=0.1$ ?SJ ま取れなかったのでとして計算$\tilde{I}_{1}(t)$$\tilde{I}_{1}(t)$$\Delta t=0.\mathrm{O}1$

Fig 2 で見ると違いがわからないが, Fig 3 では Splitting 法の誤差は相対的に Runge-Kutta 法の誤差よりもかなり小さいことがわがる.
した.

$\mathrm{t}\triangleright \mathrm{t}$’$.\cdot.\cdot’.\cdot’.\cdot.\prime^{\prime’}-\cdot’--\sim-----------------\cdot---------------------\cdot.\_{S-------}^{3-}.\cdot-\cdot----------$

.

0

.

$.\cdot.’.\cdot.\cdot.’$$\cdot$$\prime’$’.
.’$\prime’$$.2\alpha 1$’$\prime\prime\prime ’$$\prime\prime’i$$\underline{\mathrm{o}}$$.’$$.\infty\cdot\uparrow$

’

.

$i$

’
’

$j$

$\acute{\prime}$

.

$.j$

.’

$|\iota_{\mathrm{I}}.\}||\mathrm{i}|\mathrm{I}|$.

$\prime’$1$|i\mathrm{i}!.\cdot.\cdot’.\cdot’$.

$.5\triangleright’\uparrow$$\mathrm{i}_{\acute{\prime}}’\mathrm{t}|{}^{\mathrm{t}}\dot{j}$1$.\triangleright 110$

$x$

40

$\epsilon 0$$\epsilon 0$$\prime 00\mathrm{T}$Fig.2$I_{1}(t)$

の時間変化

Fig.3

$I_{1}(t)$

’20

’40

の誤差の時間変化’60.
$’\infty$

200

244
また,

と $\tilde{I}_{20}(t)=\tilde{I}_{2}(t)-I_{2}(0)\exp(-\gamma t)$ をプロットしたの力:それぞれ, Fig 4, Fig 5 である. これらの結振幅が減衰していく過程において Splitting 法の方が Runge-Kutta 法よりも正確に再現できて るこ$\tilde{I}_{2}(t)$果から,$\mathrm{A}\mathrm{a}$とがわかる.

さらに, この数値計算で非常に特徴的なことは, Splitting 法を適用した方が (4 次の)Runge-Kutta 法をを大きくとっても安定しているということである。 上の状況では (4 次の)Splitting用いた場合に比べて$\Delta t$

法では$\Delta t$$\Delta t=0.\mathrm{O}1$

を最大 $\Delta t=0.1$ くらいに大きくしても計算ができるのに対して、 (4 次) の Runge-Kutta 法で [まを大きく取っているにもかかわらず, ことしなければうまく計算できない. Splitting 法の方が$\Delta t$のような差が出るということは Splitting 法の優位性が現れているからだと考えられる.
$.\#.--$$.J’..-/^{\prime^{\wedge--\sim\sim_{\sim}}}\backslash .-\backslash --$

–

$=$

:.’.

$\backslash _{-}$$\backslash \sim_{\sim}$$\sim\sim\sim-$

$.i$

.

$\wedge\sim\sim_{\sim}-\sim-$$i’\acute{\acute{\prime}}\prime j$

.
$.i!$

$|i$

$.\cdot\ovalbox{\tt\small REJECT}$

$i$

$i$

$’\infty.,\mathrm{o}\mathrm{r}\mathrm{X}_{r}^{r}--$$\underline{u}0_{0}..\epsilon_{0}0.l0.70.\epsilon 0.50.3\mathrm{o}\mathrm{a}_{0\mathrm{a}\alpha\infty}0.\cdot\prime\prime$$\S.arrow\prime \mathrm{c}\prime\prime\sim\prime 4\cdot\prime\prime \mathrm{a}_{i}\sim^{i^{1}}’\prime\prime\prime 0^{i},,0^{\cdot}.\cdot \mathrm{a}.\acute{\alpha}.\infty\infty.’\infty’\infty’.\alpha’ \mathrm{r}.’.\infty \mathrm{a}\mathrm{o}\frac\frac{\frac{1}{}!}{}-j\frac{!||}{\overline{!}}\dot{i}|!^{i}|-|-|i^{\dot{i}}!i^{\dot{i}}|-\dot{i}i!l\acute{\acute{\prime}}\prime j\sim_{\sim}iJ_{-\backslash ----\backslash _{-\backslash \sim}}’/^{\prime^{\wedge\sim\sim_{\sim}}}\wedge-.\#--\backslash --\mathrm{T}\sim\sim-_{\wedge\sim\sim_{\sim}}-\sim-$\’i|.
$’.\alpha$$.ji^{\dot{i}}!$$\dot{_{!}|}_{,!}$$||_{_{!},!i^{\dot{i}1}}|||.\dot{.}$.
$\prime \mathrm{a}$$\infty$$’\infty \mathrm{T}$Fig.4$I_{2}(t)$I2(t) の誤差の時間変化

Fig.5

の時間変化

このように方程式 (22) に対する数値計算において, 正確性と安定性の両面から Splitting 法の有効性が確かめられた.

4Ginzburg-Landau 方程式への応用
前節と同様にして, 今度は Splitting 法をより複雑な系である複素 Ginzburg-Landau 方程式,,$\frac{\partial u}{\partial t}+p_{1}\frac{\partial^{2}u}{\partial x^{2}}+q_{1}|u|^{2}u+r_{1}u=0$に応用することについて考察しよう $[9,10]$ .

ここで$p_{1},$$q_{1},$(27)は複素数で $p_{1}=-c_{1}-\mathrm{i}c_{2},$ $q_{1}=-c_{3}-\mathrm{i}c_{4}$ ,$r_{1}$

は実数) であるとする.
ここでも, Splitting スキームを適用するために方程式 (27) を 2 つに分ける.

$r_{1}=-c_{5}-\mathrm{i}c_{6}$($c_{\mathrm{j}}$$. \frac{\partial u}{\partial t}1p_{1}\frac{\partial^{2}u}{\partial x^{2}}=0$

(28)

,

$\mathrm{i}\frac{\partial u}{\partial t}+q_{1}|u|^{2}u+r_{1}u=0$

.

(29)

前節と同様にして, 方程式 (28) で時間発展させることが $\exp[A\Delta t]$ に対応し, 方程式 (29) で時間発展させる
ことが $\exp[B\Delta t]$ に対応すると考え, 第 2 節で述べた計算手順を適用してみる. 方程式 (27) には前節で議論した (25), (26) に当たるような, 数値計算の妥当性の目安となる適当な量がない. そこで, 本節では十分に$\Delta t$

を小さく取った時の計算結果を基準にして, その基準解との差を比較の目安にする.
本節では Splitting 法の次数を大きくしていったときに誤差が小さくなっているかを確認するための計算をおこなった. 実際の数値計算の例を見てみよう. Fig 6 は初期値を非線形シュレーデインガー方程式 $(p_{1}=1.0$ ,$q_{1}=2.0,$ $r_{1}=0.0)$の2-soliton 解 $(\eta_{1}=1.0, \xi_{1}=-0.5,\eta_{2}=0.5,\xi_{2}=0.45)$ に選んで時間発展を計算したものである.

ここでは空間領域を $-25\leq x\leq 25$ で空間メツシュを $\Delta x=0.0976$ にとり, 周期境界条件としている. また,245
時間 $0\leq t\leq 25$ に対して数値計算をおこな$\Delta t=0.\mathrm{O}\mathrm{O}1$$\mathrm{A}\mathrm{a}$, 時間メッシュは$\Delta t=0.\mathrm{O}1$にとる. 基準となる解については(4 次のスキームで計算) として計算をおこなった. また (27) のパラメータは $a_{1}=0.05,$ $a_{2}=1.0$ ,-0.005, $a_{6}=0.0$ と選んだ.

$a_{3}=0.01,$ $a_{4}=2.0,$$a_{5}=$

Splitting 法の次数を 1,2,4,6,8 と変えていって, 誤差の最大値を比較したものが表 1 である. またこれをグラフにしたものが Fig 7 である. 表 1 からは, 次数を上げるに従って誤差が減少しているということがわがり, 実際に Splitting 法の計算がほぼ次数通りにうまくいっていることが確認できる.
$0.0\uparrow$o.oo’$|\mathrm{u}|$0 屋 001$’\epsilon- 05$$\frac{-33^{1}}{\in \mathrm{a}}1\mathrm{e}46$$’\cdot- 07$$1\cdot- 09$

1e-l0

Fig 6

$\mathrm{G}\mathrm{L}$13

方程式 (27) の時間発展の様子 (2 パルス解)

5

4

6

7

$\epsilon$

Fig 7Splitting 法の次数による誤差の違い

methodmesh size[1] SP 1$\Delta t=1.0\cross 10^{-}$[2] SP 2$\Delta t=1.0\cross 10^{-2}$2500$1.480\mathrm{x}10^{-4}$[3]$\mathrm{S}\mathrm{P}4$$\Delta t=1.0\cross 10^{-2}$2500$6.681\cross 10^{-7}$[4]$\mathrm{S}\mathrm{P}6$$\Delta t=1.0\cross 10^{-2}$2500$4.004\cross 10^{-10}$time step2500$\max|u_{X}-\tilde{u}_{X}|$$2.207\cross 10^{-}$[5]$\underline{\mathrm{S}\mathrm{P}8\Delta t=1.0\cross 10^{-2}25001.183\mathrm{x}10^{-10}}$

表 1Splitting 法の次数による誤差の違い (2 パルス解)

さらに, Fig 8 は初期値を $u=\sin x$ と選んで時間発展させたときの図である. ここでは, 空間の解像度の関係
から空間メッシュを $\Delta x=0.0488^{\cdot}$ と取っている. この初期値に対する数値計算においても Splitting 法の次数を上げていくと, 表 2, Fig 9 から誤差の最大値が次数に従って減少していることが確認できる.
$000$ ’

$0D\infty 1$$|\mathrm{u}|$1会6$,*\cdot 06$$\frac{-3\dot{3}}{\in \mathrm{a}}$

’“7

’e-08

$1\cdot- \mathit{0}9$$1\mathrm{r}\prime 0$$1\cdot\prime 11$Fig 8$\mathrm{G}\mathrm{L}$

方程式 (27) の時間発展の様子 ( $\sin$ 波).

5$\mathrm{S}1\emptyset*-$Fig 9Splitting 法の次数による誤差の違い$\epsilon$

246
methodth[1] SP$\mathrm{d}$mesh size$\mathrm{h}$$\mathrm{i}$time step$\max|u_{X}-\tilde{u}_{X}|$$1\Delta t=1.0\mathrm{x}10^{-}25001.9\ovalbox{\tt\small REJECT} 38\cross 10-$[2]$\mathrm{S}\mathrm{P}2$$\Delta t=1.0\cross 10^{-2}$2500$4.746\cross 10^{-5}$[3]$\mathrm{S}\mathrm{P}4$$\Delta t=1.\mathrm{O}\cross 10^{-2}$2500$2.142\mathrm{x}10^{-6}$[4]$\mathrm{S}\mathrm{P}6$$\Delta t=1.\mathrm{O}\mathrm{x}10^{-2}$2500$2.864\cross 10^{-9}$

表 2Splitting 法の次数による誤差の違い ( $\sin$ 波)

方程式 (27) についても, これらの例のように Hamiliton 系からの差があまり大きくない場合はSplitting 法が有効で, しかもその次数を上げると数値誤差 (打ち切り誤差) が小さくなっていく様子が確認で従って,$\mathrm{G}\mathrm{L}$きた.

また, この計算においても前節同様 (4 次の)Splitting 法は (4 次) の Runge-Kutta 法に比較してきく取っても安定であることが確認できた.

$\Delta t$

を大

5 まとめ
本稿では, Hamilton 系に属する非線形波動方程式 (14) に対する Symplectic 数値解法と, 方程式 (14) に散逸項や拡散項を加えた方程式, (22),(27) に対する Splitting 法の有効性について考察をおこなった. 特に,
これまで不安定になると言われていた 3 次以上の Splitting 法を方程式 (27) に適用し, パラメータが限られた範囲ではあるが, 不安定になることなくほぼ次数通りの計算精度が出ることを確認した.
本研究は次の 2 点で, 特に意義があると考えられる. まず, 今回考察をおこなった全ての系において Splittingを大きく取れることがわかった.
法あるいは Symplectic 法は (4 次の古典)Runge-Kutta 法よりも安定で$\Delta t$

Runge-Kutta 法は広く用いられているのでこれよりも安定で計算効率の良いスキームを挙げておくことは意
味のあることだと考えられる. 第 2 に, 陽的な Symplectic 数値解法と Splitting 法は形式的には全く同じなので, 本稿で挙げた例のように散逸項や拡散項が比較的小さい場合には計算手法を変えることなく , 方程式 (14)と方程式 (27) を統一的に扱えることがわかった. 第 1 節でも述べたように, 系が変わるごとに数値計算スキー

ムを変えなければならないとなると不便であるが, 本稿のように両方を同じ形式で扱えれば非常に都合がよい.
本研究では, 幾つかの例について Splitting 法の有効性を示した. しかし, 一般論として例えば, パラメー
タのどの範囲まで Splitting 法が安定で有効なのか等は未だ明かではない. このことが, 今後の課題だと考え
られる.

参考文献
[1] 大貫義郎, 吉田春夫, 力学 (現代物理学叢書), 岩波書店 (2001)[2] E. Hairer,.
Springer (1993).
$\mathrm{S}.\mathrm{P}$[3]$\mathrm{J}.\mathrm{M}$, G. Wanner, Solving Ordinary Differentid Equahons (2nd revised edition),$\mathrm{N}\emptyset \mathrm{r}\mathrm{s}\mathrm{e}\mathrm{t}\mathrm{t}$. Sanz-Serna, and$\mathrm{M}.\mathrm{P}$. Calvo,$Numer\dot{\tau}cal$

Hamiltonian Problems, Chapman and Hall (1994).

[4] H. Yoshida, Recent progress in the theory and application of symplectic integrators, Celest. Mech. 56(1993) 27-43.
[5] 吉田春夫, シンプレクティック数値解法, 数理科学 384 (1995) 37-46.

247[6] 佐々成正, 吉田脊夫, 非線形 Schr\"odinger 方程式に対する symplectic 数値解法応用数理 10 No 2$(20\mathrm{t}$119-131.
[7] H. Yoshida, Construction of higher order symplectic integrators, Phys. Lett A150 (1990) 262-268.
[8] A. Torcini, H. Frauenkron and P. Grassberger, Physica $D103$ (1997) 605.

[9] D. Goldman, and L. Sirovich, Quart. Appl. Math. 53 (1995) 315.

[10] D. Goldman, and$\mathrm{T}.\mathrm{J}$. Kaper, SIAM J. Numer. Anal. 33 (1996) 349.