数理解析研究所講究録
第 1543 巻 2007 年 245-253

245

Landau-Lifshitz 方程式に対する差分スキームについて
武蔵野大学 秋山高宏 (Takahiro Akiyama)lMusashino university
岐阜大学 石渡哲哉 (Tetsuya Ishiwata)2
Gifu university
みずほ情報総研 不破敦 (Atsushi Fuwa)

Mizuho Information&Research Institute早稲田大学 堤正義 (Masayoshi Tsutsumi)
Waseda university
概要

強磁性体中のスピンの運動を記述するランダウ. リフシッツ方和式に対する差分スキームを提案する。 この方程式は(1) 解の長さを保存する、(2) エネルギー保存系あるいは散逸系である、という性質をもつ。本講演では、提案する差分スキームがこの両方の性質を元の方程式から継承していることを示す。 またそのスキームは陰的非線形であるので、 差分解の一意可解性等についても言及する。 最後に、 厳密解に対する数値計算結果を紹介する。

1Landau-Lifshitz 方程式について強磁性体の中の電子スピンの動きを考える。歳差運動による制動を考慮した電子スピンの運動モデルとして、 Landau-Lifshitz 方程式$\frac{\partial u}{\partial t}=u\cross\Delta u-\mu ux(u\cross\Delta u)$が知られている。 ここで,
$u_{2}(x, t),u_{3}(x, t))$

た.

$\Delta=\sum_{:=1}^{N}\frac{\partial^{2}}{\partial x_{1}^{2}}$:$\Omega\subseteq \mathbb{R}^{N}$$\Omega x(0, \infty)arrow \mathbb{R}^{3}$であり、$\cross$(1)は強磁性体領域を表し、 $u=u(x, t)=(u_{1}(x, t)$ ,は, 電子スピンの向きを表すベクトルとする。 まは外積を表すものとす乱 この方程式の右辺第二項は Gilbert

制動項と呼ばれ、係数 は非負とする。 $\mu=0$ の場合は歳差運動から生じる制動を考慮していないモデルであり、 Heiscnberg 方程式と呼ばれる。
以下、 Landau-Lifshitz 方程式の解がもつ重要な性質として、 長さとエネルギーに関する性質について述べる。初期値 $v_{0}(x)$ が $|u_{0}(x)|=1$ を満たすとすると、解は $t>0$ で常に$|u(x, t)|=1$ を満たす。つまり、長さが保存する。 これは、方程式と $u(x, t)$ との内積をとる$:=\{(x, y, z):x^{2}+y^{2}+z^{2}=1, x, y, z\in \mathbb{R}\}$
ことによりわかる。以下、解は単位球面上$\mu$$S^{2}$

1 非常勤講師

2e.mail: tisimata@gifu-u.ac.jp

246がに値をとるもののみを考える。 ここで、 自明解について述べる。三つ組みの数$\alpha^{2}+\beta^{2}+\gamma^{2}=1$ を満たすとすると、 $u_{a}=(\alpha, \beta, \gamma)\in S^{2}$ は、 適切な境界条件の下で方程式 (1) の自明解となる。
次に、 以下で定義されるエネルギー $E(u(t))$ について述べよう:$\alpha,$$E(u(t))$$\beta,$$\gamma$$:=\Vert\nabla u(\cdot,t)\Vert_{L^{2}(\Omega)}^{2}$

このエネルギーに関しては、 方程式ととの内積を領域まで積分することにより、次の評価を得る:$\Delta u$$\Omega$及び時刻$E(u(t))=E(u_{0})-2 \mu\int_{0}^{t}||u(\cdot, s)x\Delta u(\cdot, s)||_{L^{2}(\Omega)}^{2}ds$$0$

から任鳶の $t>0$.

(2)

これより、 Landau-Li&hitz 方程式の解のエネルギーは時間 に関して非増加であることが$t$分かる。特に、 $\mu=0$ の場合はエネルギーが保存する。
本小論の目的は、 Ltdau-Li&hitz 方程式の差分スキームを構成することにある。 これまでにも既に、差分解の長さを保存する方法 ( の projection mcthod) などが知られている。 この方法は、差分スキームを用いて得たベクトルを、 その長さで割ることにより単位球面上に射影する方法であり、常に差分解の長さは 1 となる。つまり、差分解の大きさが発散することなく安定に数値計算が行えることがわかる。 しかし、残念ながらエネルギー構造について考慮されておらず、 例えば、 $\mu=0$ の Hciscnbcrg 方程式の場合にはエネルギーの保存性が崩れる。 また、 $\mu>0$ の場合は、 一見エネルギーが減少し解の空間的な起伏が穏やかになっていくため数値計算が安定に、 かつうまく行われているように見える。 しかし、 これは dumping が効いているので解があまり激しい挙動を示さないことに起因しているだけであり、厳密解と比較するとエネルギーの挙動などは異なることが分かる。エネルギー構造に着目した離散化としては、 降旗氏らの一連の研究がある。 (例えば、$[2, 6]$ など。) 彼らの研究により、 エネルギーの散逸性や保存性を継承する差分スキームの有効性が指摘されている。以上を受け、 本小論では Landau-Lifshitz 方程式がもつ、 解の長さ保存とエネルギー散逸性 (Heisenberg 方程式の場合は、 保存性) を継承する差分スキームを提案する。 ただし、 以下では簡単のため $\Omega=[-1,1]$ とし、 周期境界条件のみを考える。以下の結果は、高次元の場合や他の境界条件の場合にも適切な設定の元で成立することを注意しておく。 2 章では、 我々の提案する差分スキームを紹介し、 差分解の長さが保存すること、 および差分解のエネルギー評価について述べる。 ただし、提案する差分スキームは陰的かつ非線形であるので、 3 章ではその一意可解性について述べる。 4 章では、 提案した差分スキームを用いて、 厳密解に対する数値計算を行った結果を紹介する。
最後に、今後の課題等についてコメントする。
$[1|$$s$2

提案する差分スキームと差分解の性質

空間刻みおよび時間刻みは等間隔とし、 それぞれの刻み幅をで表す。 軸の番目の格子点を$m$ ステップ目の時刻を偏とし、 $(x, t)=(x_{n}, t_{m})(n=0,1,$ $\cdots N-1$$\Delta x,$$\Delta t$

$x$

$n$

$x_{n},$

3 もちろん、 差分スキ–ムにより計算するベクトルの長さが$0$にならない、 といった保証は必要である。

247

and $m=0,1,2,$ ) における差分解を我々の提案する差分スキームは以下である:$\cdots$$u_{n}^{m}=(u_{1,n}^{m}, u_{2,\mathfrak{n}}^{m}, u_{3,n}^{m})$と表記する。

$\frac{u_{n}^{m+1}-u_{n}^{m}}{\Delta t}=u_{n}^{m+1/2}x\Delta_{h}u_{n}^{m+1/2}-\mu u_{n}^{m+1/2}\cross(u_{n}^{m+1/2}x\Delta_{h}u_{n}^{m+1/2})$

.

(3)

境界条件は
$u_{0}^{m}=u_{N}^{m}$,$u_{-1}^{m}=u_{N-1}^{m}$$(m=0,1,2, \cdots)$(4)とし、 初期値は$(n=0,1,2, \cdots N)$$u_{n}^{0}=u_{0}(x_{n})$とする。 ここで、 $u_{n}^{m+1/2}=(u_{\mathfrak{n}}^{m+1}+u_{n}^{m})/2$ であり、$\Delta_{h}$(5)は次の二階差分作用素である:$\Delta_{h}u_{n}=\frac{u_{n+1}-2u_{n}+u_{n-1}}{\Delta x^{2}}$.

以下で、 この差分スキームが長さ保存を満たし、かつ、エネルギーの散逸性あるいは保存性を継承していることを示す。 ここで、「継承」 と書いたのは、 エネルギーの散逸性といっても、単にエネルギーが減少していることを示すのではなく、以下で導入する差分解に対するエネルギーが、元のエネルギー評価式 (2) の離散化に対応した式を満たすからである。
まず、 長さ保存についてであるが、 式 (3) ととの内積をとることにより$u_{n}^{m+1/2}$$|u_{\mathfrak{n}}^{m+1}|=|u_{n}^{m}|$が得られる。 これにより、 差分解の最大値ノルムの有界性が得られる。
次に、差分解のエネルギー $E_{h}(u^{m})$ を次で定義する:$E_{h}(u^{m})$$:=||D^{+}u^{m}||_{2}^{2}$.

ここで、 $D+$ は一階前進差分作用素:$D^{+}u_{n}$であり、 ノルム$||\cdot||_{2}$$:= \frac{u_{\mathfrak{n}+1}-u_{n}}{\Delta x}$

は$||v||_{2}=( \sum_{n=0}^{N-1}|v_{n}|^{2}\Delta x)^{1/2}$と定義する。 また、$D^{-}$を一階後退差分作用素:$D^{-}u_{n}$$:= \frac{u_{n}-u_{n-1}}{\Delta x}$.とすると、 部分和文公式:$\sum_{n=0}^{N-1}(D^{-}a_{n})b_{\mathfrak{n}}=-\sum_{n=0}^{N-1}a_{n}(D^{+}b_{n})$248に対して成立する。以上を用いると、式 (3) とに関して周期 $N$ をもつとの内積に、 上の部分和文公式を適用し、 時間方向に和をとることによりが$n$$a_{n},$$b_{n}$$\Delta_{h}u_{n}^{m+1/2}$.

$E_{h}(u^{M})=E_{h}(u^{0})-2 \mu\sum_{m=0}^{M-1}||u^{m+1/2}\cross\Delta_{h}u^{m+1/2}||_{2}^{2}\Delta t$を得る。 この式は、 式 (2) の離散化に相当することがわかる。 つまり、 スキーム (3) による差分解は正しくエネルギー構造を継承していることがわかる。

3$u_{n}^{m+1}$の一意可解性に関して陰的かつ非線形である。 よって、本前節で提案した差分スキーム (3) は、に対する一意可解性について説明する。
節ではまず、次の一意性を示す。
$u_{n}^{m+1}$$u_{\mathfrak{n}}^{m+1}$

Lemma 1 Suppose that there eansts solution of (3). If$\frac{\Delta t}{\Delta x^{2}}<\frac{1}{2+3\mu}$,then the solution is unique.
が共に $u_{n}^{m+1}=v_{n},$ $u_{n}^{m+1}=w_{n}$ として式 (3) を満たすとする。 ここで、 $a_{n}x\Delta_{h}a_{n}=a_{n}x(a_{n+1}+a_{n-1})/\Delta x^{2}$ であることに注意し、 また、

Sketch of proof.

$||(v+u^{m})/2||_{\infty}\leq 1$$||w-v||_{\infty}$$v_{n}$

および$w_{n}$および $||(w+u^{m})/2||_{\infty}\leq 1$ であるので、$\frac{\Delta t}{\Delta x^{2}}(2||\frac{w-v}{2}||_{\infty}\cdot||\frac{w+u^{m}}{2}||_{\infty}+2||\frac{v+u^{m}}{2}||_{\infty}\cdot||\frac{w-v}{2}||_{\infty}$$\leq$$+2 \mu||\frac{w-v}{2}||_{\infty}\cdot||\frac{w+u^{m}}{2}||_{\infty}^{2}+2\mu||\frac{v+u^{m}}{2}||_{\infty}\cdot||\frac{w-v}{2}||_{\infty}\cdot||\frac{w+u^{m}}{2}||_{\infty}$$+2 \mu||\frac{v+u^{m}}{2}||_{\infty}^{2}$$\leq$.

$|| \frac{w-v}{2}||_{\infty})$$(2+3 \mu)\frac{\Delta t}{\Delta x^{2}}||w-v||_{\infty}$を得る。 ここで、 $||a||_{\infty}= \max_{n}|a_{n}|$ である。 よって、 主張は証明された。
次に、 解の存在を示す。

Theorem A If$\frac{\Delta t}{\Delta x^{2}}<\min\{\frac{M-1}{M^{2}+\mu M^{3}},$for some $M>1$ , then the solution of (3) exists.

$\frac{1}{2M+3\mu M^{2}}\}$

(6)

249

証明の outline を示す。以下のようにして近似列を構成し、 この近似列が収束することを示す。
$\frac{v_{n}^{k+1}-u_{n}^{m}}{\Delta t}$$=$$\frac{v_{n}^{k}+u_{\mathfrak{n}}^{m}}{2}x\Delta_{h^{\frac{v_{n}^{k}+u_{n}^{m}}{2}}}$(7)$- \mu\frac{v_{n}^{k}+u_{n}^{m}}{2}x(\frac{v_{\mathfrak{n}}^{k}+u_{n}^{m}}{2}x\Delta_{h}\frac{v_{n}^{k}+u_{\mathfrak{n}}^{m}}{2})$

次の補題により、 近似列$\{v_{n}^{k}\}$の有界性を示す。

Lemma 2 Let $M$ be a positive constant larger than 1. Suppose that $|u_{n}^{m}|=1$ . If
(8)$|| \frac{v_{n}^{0}+u_{n}^{m}}{2}||_{\infty}\leq M$

and$\frac{\Delta t}{\Delta x^{2}}\leq\frac{M-1}{M^{2}+\mu M^{3}}$,(9)then we have$|| \frac{v_{n}^{k}+u_{n}^{m}}{2}||_{\infty}\leq M$

for any $k\in N$ .
Sketch of proof.
$w_{n}^{k+1}$$=$$w_{n}^{k}=(v_{n}^{k}+u_{n}^{m})/2$とおく。式 (7) より$\frac{v_{n}^{k+1}-u_{n}^{m}}{2}+u_{n}^{m}$$u_{\mathfrak{n}}^{m}+ \frac{\Delta t}{2\Delta x^{2}}(w_{n}^{k}x(w_{\mathfrak{n}+1}^{k}+w_{n-1}^{k})-\mu w_{n}^{k}\cross(w_{\mathfrak{n}}^{k}x(w_{n+1}^{k}+w_{n-1}^{k})))$

を得る。 よって、$||w^{k+1}||_{\infty} \leq 1+\frac{\Delta t}{\Delta x^{2}}(||w^{k}||_{\infty}^{2}+\mu||w^{k}||_{\infty}^{3})$である。 ここで、 $||w^{k}||_{\infty}\leq M$ であれば、 条件 (9) より$||w^{k+1}||_{\infty}$$M1+ \frac{\Delta t}{\Delta x^{2}}(M^{2}+\mu M^{3})$$\leq\leq$となる。 よって、結論を得る。
Remark: 反復の初期値としてする。

$v_{n}^{0}=u_{n}^{m}$とすれば、 条件 (8) は任意の $M>1$ について成立

250

次に、 近似列行ない、$\frac{v_{n}^{k+2}-v_{n}^{k+1}}{\Delta t}$$\{v_{n}^{k}\}$$=$が収束列であることを示す。 前の二つの証明と似たような式変形を$\frac{v_{n}^{k+1}-v_{\mathfrak{n}}^{k}}{2}x\frac{w_{n+1}^{k+1}+w_{n-1}^{k+1}}{\Delta x^{2}}+w_{n}^{k}x\frac{v_{\mathfrak{n}+1}^{k+1}-v_{n+1}^{k}+v_{n-1}^{k+1}-v_{n-1}^{k}}{2\Delta x^{2}}$$- \mu(\frac{v_{\mathfrak{n}}^{k+1}-v_{n}^{k}}{2}x(w_{n}^{k+1}x\frac{w_{\mathfrak{n}+1}^{k+1}+w_{\mathfrak{n}-1}^{k+1}}{\Delta x^{2}})+w_{n}^{k}x(\frac{v_{n}^{k+1}-v_{\mathfrak{n}}^{k}}{2}\cross$$\frac{w_{\mathfrak{n}+1}^{k+1}+w_{n-1}^{k+1}}{\Delta x^{2}})+w_{n}^{k}x(w_{n}^{k}x\frac{v_{\mathfrak{n}+1}^{k+1}-v_{n+1}^{k}+v_{n-1}^{k+1}-v_{n-1}^{k}}{2\Delta x^{2}}))$を得る。 前の補題より得られる近似列の有界性から、$||v^{k+2}-v^{k+1}||_{\infty} \leq(2M+3\mu M^{2})\frac{\Delta t}{\Delta x^{2}}||v^{k+1}-v^{k}||_{\infty}$である。条件 (6) より、$\{v_{n}^{k}\}$は極限値$v_{n}^{*}= \lim_{karrow\infty}v_{n}^{k}$を持つ。 よって、定理が証明された。

4

数値計算結果

この節では、提案した差分スキームによる数値計算例を紹介する。解析解と数値解のエネルギーの挙動を比較するため、 厳密解に対する数値計算を行った。厳密解は、 $\mu=0$ の場合は既に [4] により得られている解を用い、 $\mu>0$ の場合については以下で紹介する厳密解を用いた。

4.1

厳密解であるとする。 ただし、とし、Landau-Lifshitz 方程式 (1) の非自明解として以下の解がある:$\alpha,$$k\in \mathbb{R}$$u(x,t)$$\Omega=\mathbb{R}/2\pi k$$=$$\alpha\neq l\pi(l\in \mathbb{Z})$とする。 このとき、$(u_{1}(x,t),u_{2}(x, t),u_{3}(x,t))$$( \frac{\sin\alpha\cos[k\cdot x-\phi(x,t;\alpha,k,\mu)]}{d(t;\alpha,k,\mu)},$ $\frac{\sin\alpha\sin[k\cdot x-\phi(x,t;\alpha,k,\mu)]}{d(t;\alpha,k,\mu)}$$\frac{e^{k^{2}\mu l}\cos\alpha}{d(t;\alpha,k,\mu)})$.

ここで、 $d(t;\alpha,k,\mu)$ および $\phi(x,t;\alpha,k,\mu)$ はそれぞれcos2$d(t;\alpha,k,\mu)=\sqrt{\sin^{2}\alpha+e^{2k^{2}\mu t}}$$\alpha$251$\phi(x, t;\alpha, k,\mu)=\frac{1}{\mu}$log$( \frac{d(t;\alpha,k,\mu)+e^{k^{2}\mu t}\cos\alpha}{1+\cos\alpha}I$である。 この解は、 容易に高次元の場合に拡張できることを注意しておく。 4
ここで、 $\muarrow 0+$ とすると、$u_{1}(x, t)$$=$$u_{2}(x, t)$$=$$u_{3}(x, t)$$=$sin cos ( $k\cdot x-(|k|^{2}$ cos ) )sin sin ( $k\cdot x-(|k|^{2}$ cos ) )cos$\alpha$$\alpha$$\alpha$$\alpha$$t$$t$$\alpha$をえる。 このとき関数 $u(x, t)=(u_{1}(x,t),$ $u_{2}(x, t),$ $u_{3}(x,t))$ は Hciscnberg 方程式の厳密解になっている。 この解は既に [4] において得られている厳密解である。

4.2

数値実験結果図 1 と図 2 の数値実験では、 $\Delta x=1/30,$ $\Delta t=1/20000$ として、初期値として前節の厳密解において $t=0$ を代入した値を使っている。
図 1 は、 $\mu=0$ 即ち Hcnscnbcrg 方程式に対する数値実験結果である。 (a) は数値解の第一成分、 (b) は厳密解の第一成分の時間発展をそれぞれ表示している。 (c) はエネルギーの挙動である。 エネルギーがよく保存されており、 解の形状変化も相違ないことがわかる。
次に、 図 2 では、 $\mu=0.1$ の場合の数値実験結果を示す。先ほどと同様に、 (a), (b) はそれぞれ数値解および厳密解の第一成分の時間発展を表示している。数値解は厳密解の挙動をよく近似していることがわかる。 (c) は、 数値解と厳密解のエネルギーの挙動の比較で.

ある。実線は数値解を、点線は厳密解を表している。数値解は、厳密解のエネルギー散逸性をよく近似していることがわかる。
厳密解に対する数値実験のみ提示したが、数値解は以上のようによく解析解の性質を引き継いだ時間発展をしていることがわかる。

5終りに今回、方程式がもつ性質の中で長さ保存およびエネルギー構造に着目し、 それらを継承する差分スキームを提案し、厳密解に対する数値計算によりその有効性を見た。 まず思いつく課題として、 この離散化に対する誤差評価がある。 これについては最近 Sobolev ノルムによる誤差評価を終えたので、近日中に報告できれば、 と思っている。次に、今回提案したスキームは陰的かつ非線形であるので、近似列を構成する際の反復法の改良が考えられる。高速化や線形多段化などを検討する余地が大いにあるだろう。
4 この解は新たに我々が構成した解であるが、 もし既に結果として公表されていることをご存じの方はこ一報頂きたい。なお、 [3] において $Lak_{8}hmanan$ と $N\sim ura$ は Landau-Li&hitz 方程式に対する厳密解の構成法を提案しているが、彼らの主張は一般には正しくない。 ([5] を参照。) 彼らの方法では、 ここで紹介した厳密解 (10) は得られないことを注意しておく。

252(b) 厳密解の第一成分の時間発展(a) 数値解の第一成分の時間発展$rightarrow$$m$$sn$$rightarrow$$rightarrow$,.

.

.

.

(c) 数値解のエネルギーの挙動

図 1: $\mu=0$ の場合の数値計算結果例

謝辞
本小論は第二著者が京都大学数理解析研究所研究集会「波動現象の数理と応用」 にて講演した内容に基づいている。講演の機会を与えてくださった田中光宏先生には改めてお礼申し上げます。

参考文献[1] WEINAN $E,$ $XlAO$-PING WANG, Numerv cal $rn,p,th$ ) for the Landau-Lifshitz
tion, SIAM J. NUMER. ANAL. Vol. 38, No. 5, pp. 1647-1665.
$fs$$c^{J}qmz-$

That inherit Energy[2] D. FURIHATA, Finite Difference Schemes forConservation or Dissipation Property, J. Comput. Phys. 156 (1999), 181-205$\frac{\partial u}{\partial t}=(\frac{\partial}{\partial x})^{\alpha}\frac{\delta G}{\delta u}$[3] M.LAKSHMANAN, K.NAKAMURA, Landau-Lifshitz Equation of Ferromagnetism:
Exact ffiatment of the Gilbert Damping, Physical Rcvicw Lcttcrs 53, 2497-2499(1984)

253(a) 数値解の第一成分の時間発展(b) 厳密解の第一成分の時間発展(c) 数値解 (実線) と厳密解 (点線)のエネルギーの挙動の比較

図 2: $\mu=0.1$ の場合の数値計算結果例[4] M.LAKSHMANAN, TH.W.RUIJGROK AND C.J.THOMPSON, On the dynamicals of acontinuum $sp$ in system, Physca A 84 (1976), 577-590.
[5] E. MAGYARI, H. THOMAS, R. WEBER, Comment on “Landau-Lifshitz Equation ofFerromagn etism: Exact Truatment of the Cilbe Damping“, Physical Review Letters56, 1756 (1986)$rt$[6] T. Matsuo and D. Furihata, Dissipativ , or $Conse,rvati\uparrow$) Finite-Differen ce Schemes
for Complex- Valued Nonlinear Partial Differential Equations, J. Comput. Phys. 171(2001),.
$J$$425\triangleleft 47$$e$