令和 2 年度
卒業論文
時間発展方程式の解法に用いる
アフィン数値積分法に関する研究

川畑唯一(学籍番号：15232025)

九州工業大学

情報工学部

電子情報工学科
指導教員：小田部荘司

令和 3 年 2 月 16 日

教授

目次
第 1 章 序論 .............................................................................................................................2
1.1 はじめに ........................................................................................................................2
1.2 第一種超伝導体と第二種超伝導体 ...............................................................................2
1.3 (Time-Dependent) Ginzburg-Landau 方程式.............................................................3
1.3.1 Ginzburg-Landau(GL) 方程式 ................................................................................3
1.3.2 Time-Dependent G-L (TDGL) 方程式...................................................................8
1.3.3

TDGL 方程式の簡易化 ........................................................................................9

1.4 アフィン数値積分法(Affine Integrator(AFI)) ............................................................13
1.4.1

AFI のイメージ ..................................................................................................13

1.4.2 1 次元空間において AFI 法の時間発展方程式の導出 ..........................................14
1.4.3 2 次元空間において AFI 法の時間発展方程式の導出 ..........................................19
1.5 Lie-Trotter-Suzuki(LTS)分解 ....................................................................................21
1.6 研究の目的 .................................................................................................................23
第 2 章 実装方法と計算条件 ..................................................................................................24
2.1 実装方法.....................................................................................................................24
2.2 1 次元拡散方程式の誤差評価 .....................................................................................24
2.3 2 次元拡散方程式の誤差評価 .....................................................................................26
2.4 1 次元量子拡散方程式の誤差評価(ディリクレ境界条件) ...........................................27
2.5 1 次元量子拡散方程式の誤差評価(周期的境界条件) ..................................................29
2.6 1 次元 TDGL 方程式の誤差評価................................................................................31
2.7 1 次元 TDGL 方程式の安定性 ...................................................................................32
第 3 章 実装結果および考察 ..................................................................................................33
3.1 1 次元拡散方程式の誤差評価 .....................................................................................33
3.2 2 次元拡散方程式の誤差評価 .....................................................................................34
3.3 1 次元量子拡散方程式の誤差評価(ディリクレ境界条件) ...........................................35
3.4 1 次元量子拡散方程式の誤差評価(周期的境界条件) ..................................................36
3.5 1 次元 TDGL 方程式の誤差評価................................................................................37
3.6 1 次元 TDGL 方程式の安定性 ...................................................................................41
第 5 章 結論 ...........................................................................................................................43
謝辞 ........................................................................................................................................45
研究業績.................................................................................................................................46
参考文献.................................................................................................................................47

1

第 1 章 序論
1.1 はじめに1911 年に Kamerlingh Onnes が液体ヘリウムを用いて水銀を冷却することにより、電気抵抗がゼロになることを発見した。1933 年には、W. Meissner と R. Ochsenfeld によって完
全反磁性(マイスナー効果)という超伝導体内の磁束密度はゼロに保持されるという性質を発見した。このような「電気抵抗がゼロ」や「完全反磁性(マイスナー効果)」という現象を超伝導現象といい、超伝導現象が起きている状態を超伝導状態、そうでない状態を常伝導状態と呼ぶ。超伝導現象を持つ物質を超伝導体と呼ぶ。この発見の後、多くの物質からも超伝導現象が確認された。超伝導体は電気抵抗ゼロという現象から多くの工学的応用が期待された。しかし、その性質は温度や磁界、電流の変化によって、失われてしまう。常伝導状態から超伝導状態に転移する温度や磁界、電流をそれぞれ臨界温度 𝑇 、臨海磁界 𝐻 、臨界電流 𝐽 と呼ぶ。
1957 年に J. Bardeen、L. N. Cooper、J. R. Schrieffer らによって BCS 理論が提唱され、超伝導現象の基本的な発現機構が明らかになりました。このとき BCS 理論において 𝑇 が30 T を超えないと予想されていた。しかし、1986 年に J. G. Bednorz や K. A. Muller によって𝑇 が 35 K となる La-Ba-Cu-O 系の超伝導体が発見された。この発見以降、世界各国で𝑇 の高い超伝導体の探索が行われ、1987 年に窒素の沸点 (77 K)より高い𝑇 をもつ超伝導体が発見された。このことで入手が困難なヘリウムではなく窒素を冷媒として用いることができるようになったため、冷却コストの低減により様々な機器への応用が期待されている。
しかしこれらの高温超伝導体にも実用化に向けて様々な課題が残されているため現在でも研究が続けられている。

1.2 第一種超伝導体と第二種超伝導体
超伝導体は電気抵抗ゼロ、マイスナー効果の 2 つの性質を持っている。通常はマイスナ
ー効果によって超伝導内部の磁束密度はゼロに保たれているが、外部磁場𝐻 を高くしていき、臨界磁場𝐻 を超えると、マイスナー効果が失われ常伝導状態となる。超伝導体は超伝導状態から常伝導状態へと遷移する状態の違いによって第一種超伝導体、第二種超伝導体の 2種類に分類できる。
第一種超伝導体も第二種超伝導体も𝐻 が低い領域ではマイスナー効果を示す。第一種超伝導体は𝐻 以上の𝐻 をかけるとマイスナー効果を完全に失い、常伝導状態へと遷移する。
対して第二種超伝導体は下部臨界磁界𝐻 以上の𝐻 をかけたとき、超伝導体が破壊され磁束が一部侵入する。しかし、超伝導状態は完全には破壊されず、超伝導状態と常伝導状態が混
2

在した状態となる。超伝導状態と常伝導状態が混在している、𝐻 が𝐻 以上で上部臨界磁場
𝐻 未満における状態を混合状態と呼ぶ。𝐻 を高くすると常伝導状態の範囲が増え、入り込
む磁束線の数が増える。𝐻 が𝐻 以上となると超伝導状態は完全に破壊され、常伝導状態に遷移する。このように、第一種超伝導体と第二種超伝導体は混合状態の有無で分類されている。

1.3 (Time-Dependent) Ginzburg-Landau 方程式
ここでは、超伝導体の数値解析に用いた時間依存 Ginzburg-Landau(Time-Dependent
Ginzburg-Landau)方程式について説明する。

1.3.1 Ginzburg-Landau(GL) 方程式
Ginzburg-Landau(G-L)理論は、1950 年に V. L. Ginzburg と L. D. Landau によって提唱された超伝導現象を説明する現象論である[8]。G-L 理論は磁界が存在するときの超伝導が共存する場合の相転移を取り扱ったもので、特に第 2 種超伝導体の磁気特性を良く記述す
ることが知られている。以下では、G-L 理論を成り立たせる仮定とそれによって求められる方程式を説明する。
G-L 理論では、まず超伝導状態の秩序の程度を表す量として複素数であるオーダパラメ
ータ𝛹 = |𝛹|exp(i𝜑)を定義する。そして𝛹は、次の関係を満たすと仮定する。
(1.1)|𝛹| ∝ 𝑛ここで、𝑛 は超伝導電子密度である。超伝導状態の自由エネルギー𝐸 は𝑛 に依存しているから(1.1)式より、|𝛹| の関数である。よって、𝐸 は|𝛹| のべき展開で表すことが出来る。
𝐸 = 𝐸 + 𝛼|𝛹| +

𝛽
|𝛹|
2

(1.2)

ここで、𝐸 は常伝導状態の自由エネルギー、𝛼および𝛽はそれぞれべき展開した際の 1 次と 2 次の展開係数であり、𝑇 < 𝑇 では𝛼 < 0, 𝛽 > 0である。2 次の項までの展開であるのは|𝛹| が転移点近傍(𝑇 ≅ 𝑇 )において十分に小さいからである。
次に、磁界の存在が𝛹の空間的変化に寄与することを考慮し、(1.2)式に磁界のエネルギー
密度と運動エネルギー密度を加算する。
𝐸 = 𝐸 + 𝛼|𝛹| +

𝛽
1
1
|𝛹| +
(𝛁 × 𝑨) +|(−iℏ𝛁 + 𝑒 ∗ 𝑨)𝛹|
2
2𝜇
2𝑚∗

(1.3)

ここで、𝜇 は真空中の透磁率、𝑨はベクトルポテンシャル、𝑚∗は超伝導電子の質量、ℏはディラック定数、𝑒 ∗は超伝導電子の電荷量、iは虚数単位である。
𝛹と𝑨は(1.3)式の𝐸 が最小なるように決定される。よって𝛹の共役複素数𝛹 ∗と𝑨について
変分法を適用する。
3

δ𝐸
𝜕𝐸
∂𝐸
=
− 𝛁∙
=0
∗
∗
𝛿𝛹
∂𝛹
𝜕𝜵𝛹 ∗

(1.4)

δ𝐸
𝜕𝐸
∂𝐸
=
− 𝛁∙
=0
𝛿𝑨
∂𝑨
𝜕𝜵𝑨

(1.5)

まず、(1.4)式について解く。𝐸 を𝛹 ∗で偏微分すると以下のようになる。
𝜕𝐸𝜕𝛽1(𝛁 × 𝑨)=𝐸 + 𝛼|𝛹| + |𝛹| +
∗
∗
𝜕𝛹
𝜕𝛹
2
2𝜇
+
=

1
|(−𝑖ℏ𝛁 + 𝑒 ∗ 𝑨)𝛹|
2𝑚∗

𝜕
𝛽
1
(𝛁 × 𝑨)𝐸 + 𝛼|𝛹| + |𝛹| +
∗
𝜕𝛹
2
2𝜇
+

1
{ℏ 𝛁𝛹 ∙ 𝛁𝛹 ∗ + 𝑒 ∗ 𝑨𝟐 𝛹𝛹 ∗ + 𝑖ℏ𝑒 ∗𝑨
2𝑚∗

(1.6)

∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}= 0 + 𝛼𝛹 + 𝛽𝛹|𝛹| + 0+

1{0 + 𝑒 ∗ 𝑨𝟐 𝛹 + iℏ𝑒 ∗ 𝑨 ∙ (0 − 𝛁𝛹)}2𝑚∗

= 𝛼𝛹 + 𝛽𝛹|𝛹| +1(𝑒 ∗ 𝑨𝟐 𝛹 − iℏ𝑒 ∗𝑨 ∙ 𝛁𝛹)
2𝑚∗

続いて、𝐸 を𝛁𝛹 ∗で偏微分すると以下のようになる。
𝛁∙

𝜕𝐸
𝜕
𝛽
1
(𝛁 × 𝑨)=𝛁∙
𝐸 + 𝛼|𝛹| + |𝛹| +
∗
∗
𝜕𝛁𝛹
𝜕𝛁𝛹
2
2𝜇
+

1
{ℏ 𝛁𝛹 ∙ 𝛁𝛹 ∗ + 𝑒 ∗ 𝑨𝟐 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝑨
2𝑚∗∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}1{ℏ 𝛁Ψ + 0 + iℏ𝑒 ∗ 𝑨 ∙ (𝛹 − 0)}= 0+0+0+0+
2𝑚∗
=𝛁∙
=

(1.7)

1
(ℏ 𝛁𝛹 + 0 + iℏ𝑒 ∗ 𝑨𝛹)2𝑚∗1ℏ 𝛁 𝟐 𝛹 + iℏ𝑒 ∗ 𝛁(𝑨 ∙ 𝛹)
2𝑚∗

ここで、𝛁(𝑨 ∙ 𝛹)に対して積の微分公式を用いると以下のようになる。
𝛁∙𝜕𝐹1(ℏ 𝛁𝟐 𝛹 + iℏ𝑒 ∗ (𝛁 ∙ 𝑨)𝛹 + iℏ𝑒 ∗𝑨 ∙ 𝛁𝛹)
=
∗
𝜕𝛁𝛹
2𝑚∗
4

(1.8)

クーロンゲージ𝛁 ∙ 𝑨 = 0より以下のようになる。
𝛁∙𝜕𝐹1(ℏ 𝛁 𝟐 𝛹 + iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹)=𝜕𝛁𝛹 ∗ 2𝑚∗

(1.9)(1.6)式と(1.9)式から(1.4)式は以下のようになる。
𝛼𝛹 + 𝛽𝛹|𝛹| +1(𝑒 ∗ 𝑨𝟐 𝛹 − iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹)2𝑚∗−1(ℏ 𝛁𝟐 𝛹 + iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹) = 0
2𝑚∗1(−ℏ 𝛁𝟐 𝛹 − 2iℏ𝑒 ∗ 𝑨 ∙ 𝛁𝛹 + +𝑒 ∗ 𝑨𝟐 𝛹) = 0
𝛼𝛹 + 𝛽𝛹|𝛹| +2𝑚∗
𝛼𝛹 + 𝛽𝛹|𝛹| +(1.10)1(−ℏ 𝛁𝟐 − 2iℏ𝑒 ∗ 𝑨 ∙ 𝛁 + +𝑒 ∗ 𝑨𝟐 )𝛹 = 02𝑚∗

因数分解すると最終的に次のようになる。
𝛼𝛹 + 𝛽𝛹|𝛹| +1(−iℏ𝛁 + 𝑒 ∗ 𝑨) 𝛹 = 02𝑚∗(1.11)

次に、(1.5)式について解く。𝐸 を𝑨で偏微分する部分を考える。𝑨はベクトルであるので、スカラー𝐸 を𝑨で偏微分したものを次のように要素ごとに表現できる。
𝜕𝐸
𝜕𝐸
𝜕𝐸
𝜕𝐸
=
𝒊𝒙 +
𝒊𝒚 +
𝒊
𝜕𝑨 𝜕𝐴
𝜕𝐴
𝜕𝐴 𝒛

(1.12)𝐴 , 𝐴 , 𝐴 はそれぞれ𝑨の𝑥, 𝑦, 𝑧成分である。ここで(1.12)式の𝑥成分に注目すると以下のよ
うになる。
𝜕𝐸
𝜕
=
𝜕𝐴
𝜕𝐴

𝐸 + 𝛼|𝛹| +
1
+
2𝑚∗

=

𝜕
𝜕𝐴

0+0+0+
+

=

𝛽
1
|𝛹| +
(𝛁 × 𝑨)22𝜇+ℏ 𝛁𝛹 ∙ 𝛁𝛹 ∗+𝑒 ∗ 𝑨𝟐 𝛹𝛹 ∗∗+iℏ𝑒 𝑨 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗𝛁𝛹)
(1.13)

1(𝛁 × 𝑨)
2𝜇

0
1
∗
𝟐
+𝑒
𝑨
𝛹𝛹 ∗
2𝑚∗ +iℏ𝑒 ∗𝑨 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗𝛁𝛹)1 𝜕
1 𝜕
(𝛁 × 𝑨) +2𝜇 𝜕𝐴2𝑚∗ 𝜕𝐴5𝑒 ∗ 𝑨𝟐 𝛹𝛹 ∗+iℏ𝑒 𝑨 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)∗⎧+ 𝜕𝐴 − 𝜕𝐴
𝜕𝑦
𝜕𝑧
⎪
⎪
1 𝜕
𝜕𝐴
𝜕𝐴
=
+
−
2𝜇 𝜕𝐴 ⎨
𝜕𝑧
𝜕𝑥
⎪ 𝜕𝐴
𝜕𝐴
⎪+
−
𝜕𝑥
𝜕𝑦
⎩

⎫
⎪
⎪
⎬
⎪
⎪
⎭

+𝑒 ∗ 𝐴 + 𝐴 + 𝐴 𝛹𝛹 ∗
⎧
⎫
⎪
1 𝜕 ⎪
+𝐴 ∙ 𝒊𝒙+∗∗∗∗
2𝑚 𝜕𝐴 ⎨+iℏ𝑒 +𝐴 ∙ 𝒊𝒚 ∙ (𝛹𝛁𝛹 − 𝛹 𝛁𝛹)⎬⎪⎪+𝐴 ∙ 𝒊𝒛⎩⎭+𝑒 ∗ (2𝐴 + 0 + 0)𝛹𝛹 ∗
⎧
⎫
⎪
1
1 ⎪
+1 ∙ 𝒊𝒙(0 + 0 + 0) +=∗ +0 ∙ 𝒊∗∗∗
2𝜇2𝑚 ⎨+iℏ𝑒𝒚 ∙ (𝛹𝛁𝛹 − 𝛹 𝛁𝛹)⎬
⎪
⎪
+0
∙
𝒊
𝒛
⎩
⎭
=

1
{2𝑒 ∗ 𝐴 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒙 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}2𝑚∗

𝑦成分、𝑧成分も(1.13)式と同様に変形できる。
𝜕𝐸
1=2𝑒 ∗ 𝐴 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒚 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)𝜕𝐴2𝑚∗(1.14)

𝜕𝐸1{2𝑒 ∗ 𝐴 ΨΨ∗ + iℏ𝑒 ∗ 𝒊𝒛 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}
=
𝜕𝐴
2𝑚∗

(1.15)(1.13)式、(1.14)式、(1.15)式から得られる,,をすべて(1.12)式に代入すると以下のようになる。
+{2𝑒 ∗ 𝐴 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒙 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}𝒊𝒙
𝜕𝐸1+ 2𝑒 ∗ 𝐴 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒚 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹) 𝒊𝒚=
𝜕𝑨 2𝑚∗+{2𝑒 ∗ 𝐴 𝛹𝛹 ∗ + iℏ𝑒 ∗ 𝒊𝒛 ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}𝒊𝒛(1.16)

成分をまとめると次のようになる。
𝜕𝐸1{2𝑒 ∗ 𝑨𝛹𝛹 ∗ + iℏ𝑒 ∗ ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗ 𝛁𝛹)}=𝜕𝑨 2𝑚∗(1.17)

続いて、𝐸 を𝛁𝑨で偏微分したものも要素ごとに表現できる。
𝛁∙

𝜕𝐸
𝜕𝐸
𝜕𝐸
𝜕𝐸
=𝛁∙
𝒊𝒙 + 𝛁 ∙𝒊𝒚 + 𝛁 ∙
𝒊
𝜕𝛁𝑨
𝜕𝛁𝐴
𝜕𝛁𝐴
𝜕𝛁𝐴 𝒛

ここで、(1.18)式の𝑥成分に注目すると以下のように変形できる。

6(1.18)

𝛁∙

𝜕𝐸
𝜕 𝜕𝐸
𝜕 𝜕𝐸
𝜕 𝜕𝐸
=
∙
+
∙
+ ∙
𝜕𝛁𝐴
𝜕𝑥 𝜕 ∂𝐴𝜕𝑦 𝜕 ∂𝐴𝜕𝑧 𝜕 ∂𝐴∂𝑥∂𝑦∂𝑧(1.19)

ここで、𝐸 のうち微分された𝐴 を含むのは、ローテーションを持つ項すなわち(1.3)式でいうところの右辺第 4 項だけであるので、それ以外の項は 0 となり、(1.19)式は以下のようになる。
⎧+ 𝜕𝐴 − 𝜕𝐴
𝜕𝑦
𝜕𝑧
⎪
⎪
𝜕𝐹
1 𝜕
𝜕
𝜕𝐴
𝜕𝐴
𝛁∙
=
∙
+
−
𝜕𝐴
𝜕𝛁𝐴
2𝜇 𝜕𝑥 𝜕
⎨
𝜕𝑧
𝜕𝑥
𝜕𝑥 ⎪
⎪+ 𝜕𝐴 − 𝜕𝐴
𝜕𝑥
𝜕𝑦
⎩

,

,

⎧+ 𝜕𝐴 − 𝜕𝐴
𝜕𝑦
𝜕𝑧
⎪
⎪
1 𝜕
𝜕
𝜕𝐴
𝜕𝐴
+
∙
+
−
𝜕𝐴
2𝜇 𝜕𝑦 𝜕
⎨
𝜕𝑧
𝜕𝑥
𝜕𝑦 ⎪
⎪+ 𝜕𝐴 − 𝜕𝐴
𝜕𝑥
𝜕𝑦
⎩

⎫
⎪
⎪

⎧+ 𝜕𝐴 − 𝜕𝐴
𝜕𝑦
𝜕𝑧
⎪
⎪
1 𝜕
𝜕
𝜕𝐴
𝜕𝐴
+
∙
+
−
𝜕𝐴
2𝜇 𝜕𝑧 𝜕
⎨
𝜕𝑧
𝜕𝑥
𝜕𝑧 ⎪
⎪+ 𝜕𝐴 − 𝜕𝐴
𝜕𝑥
𝜕𝑦
⎩

⎫
⎪
⎪

⎫
⎪
⎪
⎬
⎪
⎪
⎭

(1.20)

⎬
⎪
⎪
⎭

⎬
⎪
⎪
⎭

による偏微分を行うと次のようになる。
+0+0𝜕𝐸1 𝜕 +01 𝜕 ⎧𝛁∙=∙ +0 +
∙
𝜕𝐴
𝜕𝐴
𝜕𝛁𝐴
2𝜇 𝜕𝑥
2𝜇 𝜕𝑦 ⎨+2
−
+0
𝜕𝑥
𝜕𝑦
⎩
+0
1 𝜕
𝜕𝐴
𝜕𝐴
+
∙ +
(+1)
−
2𝜇 𝜕𝑧
𝜕𝑧
𝜕𝑥
+0
=−

1 𝜕 𝜕𝐴𝜕𝐴−
𝜇 𝜕𝑦 𝜕𝑥
𝜕𝑦
=−

−

𝜕 𝜕𝐴
𝜕𝐴
−
𝜕𝑧 𝜕𝑧
𝜕𝑥

1 𝜕
𝜕
(rot 𝑨| ) −rot 𝑨|
𝜇 𝜕𝑦
𝜕𝑧
7

⎫
(−1)⎬
⎭

(1.21)

=−rot(rot 𝑨) |𝜇

𝑦成分、𝑧成分も(1.21)式と同様に変形できる。
𝛁∙rot(rot𝑨)|
𝜕𝐸
=−
𝜕𝛁𝐴
𝜇

(1.22)

𝛁∙𝜕𝐸
rot(rot𝑨)|
=−
𝜕𝛁𝐴
𝜇

(1.23)(1.21)式、(1.22)式、(1.23)式から得られる𝛁 ∙ 𝛁,𝛁∙ 𝛁, 𝛁 ∙ 𝛁 をすべて(1.18)式に代入すると以下のようになる。
𝛁∙rot(rot 𝑨) |𝜕𝐸
rot(rot 𝑨) |rot(rot 𝑨) |
=−
𝒊𝒙 −
𝒊𝒚 −
𝒊𝒛
𝜕𝛁𝑨
𝜇
𝜇
𝜇
rot(rot 𝑨)=−
𝜇(1.24)(1.17)式と(1.24)式から(1.5)式は次のようになる。
1rot(rot 𝑨){2𝑒 ∗ 𝑨𝛹𝛹 ∗ + iℏ𝑒 ∗ ∙ (𝛹𝛁𝛹 ∗ − 𝛹 ∗𝛁𝛹)} +
=0
∗
2𝑚
𝜇
1
iℏ𝑒 ∗ ∗𝑒 ∗𝟐∗)
(𝛹|𝛹|𝟐 𝑨
𝛁×𝛁×𝑨 =
𝛁𝛹
−
𝛹𝛁𝛹
−
𝜇
2𝑚∗
𝑚∗

(1.25)

よって(1.4)式と(1.5)式は、次の 2 式になる。
1(−iℏ𝛁 + 𝑒 ∗ 𝑨) 𝛹 + α𝛹 + 𝛽|𝛹| 𝛹 = 02𝑚∗(1.11)再掲

1iℏ𝑒 ∗ ∗𝑒 ∗𝟐∗)
(𝛹|𝛹|𝟐 𝑨
𝛁×𝛁×𝑨=
𝛁𝛹
−
𝛹𝛁𝛹
−
𝜇
2𝑚∗
𝑚∗

(1.25)再掲

これら(1.11)式と(1.25)式を、G-L 方程式と呼ぶ。

1.3.2 Time-Dependent G-L (TDGL) 方程式
Time-Dependent G-L(TDGL)方程式は、G-L方程式に時間依存性を付与したものである。
(1.4)及び(1.5)式に対して時間発展する場合を考えると、以下の2式となる。
δ𝐸
𝜕𝐸
𝜕𝐸
𝜕𝛹
=
− 𝛁 ∙
= −γ
∗
∗
∗
𝛿𝛹
𝜕𝛹
𝜕∇𝛹
𝜕𝑡

(1.26)

δ𝐸
𝜕𝐸
𝜕𝐸
𝜕𝑨
=
− 𝛁 ∙
= −ν
𝛿𝑨
𝜕𝑨
𝜕𝛁𝑨
𝜕𝑡

(1.27)

ここで、γとνはそれぞれ𝛹と𝑨の時定数である。さらにゲージ変換を与えると次のように
8

なる。
ℏ

∂𝛹
𝜕
→ ℏ + i𝑒 ∗ 𝑉 𝛹𝜕𝑡𝜕𝑡(1.28)

∂𝑨 𝜕𝑨
→
+ 𝛁𝑉
𝜕𝑡
𝜕𝑡

(1.29)

ここで、𝑉はスカラーポテンシャルである。(1.28)、(1.29)式を(1.26)、(1.27)式それぞれに代入すると最終的に次のようになる。
𝜕𝛹
1(−iℏ𝛁 + 𝑒 ∗ 𝑨) 𝛹 + α𝛹 + 𝛽|𝛹| 𝛹 = 0+ i𝑒 ∗ 𝑉𝛹 +𝜕𝑡2𝑚∗

(1.30)

𝜕𝑨1iℏ𝑒 ∗ ∗𝑒 ∗𝟐∗)
(𝛹|𝛹|𝟐 𝑨 = 0+ 𝛁𝑉 + 𝛁 × 𝛁 × 𝑨 +
𝛁𝛹
−
𝛹𝛁𝛹
−
𝜕𝑡
𝜇
2𝑚∗
𝑚∗

(1.31)𝛾𝜈

1.3.3

TDGL 方程式の簡易化

本研究ではTDGL方程式を数値解析するが、(1.30)、(1.31)式をそのまま解くのは困難であるため、細線近似と規格化による2つの簡易化を行う。
本研究の数値解析では、𝑨は時間に対して一定に与えるから(1.31)式の左辺第1項の時間偏微分が0となる。次に(1.31)、(1.32)に対して規格化を行う。超伝導体のコヒーレンス⾧𝜉と磁界侵入⾧𝜆は熱力学的臨界磁界を𝐻 として以下の2式のように表す。
𝜉=𝜆=ℏ
2𝑚∗ |𝛼|𝑒 ∗𝜇 𝐻𝑚∗|𝛼|(1.32)(1.33)

そして、以下の式に記す規格化を行う。
𝜉𝛁 → 𝛁

(1.34)|𝛼|𝑡→𝑡
𝛾(1.35)

𝑒 ∗𝛾
𝑉→𝑉
|𝛼|

(1.36)

𝜆
√2𝜇 𝐻
9

𝑨→𝑨

(1.37)

𝛽|𝛼|(1.38)

𝛹→𝛹

(1.30)式の規格化を行う。(1.30)式の左辺第1項は、次のようになる。

→𝛾

𝜕
𝛾
𝜕
𝑡
|𝛼|

𝛾

𝜕𝛹
+ i𝑒 ∗ 𝑉𝛹
𝜕𝑡

|𝛼|
𝛽

𝛹 + i𝑒 ∗

𝛾

|𝛼|
𝛽

=𝛾

|𝛼|
1
|𝛼|
𝛾
𝛽

= |𝛼|

|𝛼|
𝑉
𝑒 ∗𝛾

|𝛼| |𝛼|
𝜕𝛹
+i
𝛾
𝛾 𝛽
𝜕
𝑡
|𝛼|

|𝛼|
𝛽

|𝛼|
𝛽

𝛹

(1.39)

𝑉𝛹

|𝛼|
𝜕𝛹 1
+ |𝛼|
𝜕𝑡 𝛾
𝛽

i𝑉𝛹

𝜕𝛹
+ i𝑉𝛹
𝜕𝑡

(1.30)式の左辺第2項は、次のようになる。
1(−iℏ𝛁 + 𝑒 ∗ 𝑨) 𝛹
2𝑚∗
→

1
1
√2𝜇 𝐻
−iℏ 𝛁 − 𝑒 ∗
𝑨
2𝑚∗
𝜉
𝜆

= −iℏ

=

1
√2𝑚∗ 𝜉

|𝛼|
𝛽

= |𝛼|

−iℏ
|𝛼|
𝛽

𝛁 − 𝑒∗√2𝜇 𝐻√2𝑚∗ 𝜆2𝑚∗ |𝛼|√2𝑚∗ℏ𝑨𝛁 − 𝑒∗

|𝛼|
𝛽

𝛹

|𝛼|
𝛽

𝛹

√2𝜇 𝐻

𝑚∗ |𝛼|√2𝑚∗𝑒 ∗𝜇 𝐻

(1.40)

𝑨𝛹

(−i𝛁 − 𝑨) 𝛹

(1.30)式の左辺第3項は、次のようになる。
𝛼𝛹 → 𝛼

|𝛼|
𝛽

10

𝛹

(1.41)

(1.30)式の左辺第4項は、次のようになる。

𝛽|𝛹| 𝛹 → 𝛽

= |𝛼|

|𝛼|
𝛽

𝛹

|𝛼|
𝛽

|𝛼|
𝛽

|𝛹| 𝛹

𝛹
(1.42)(1.39)式から(1.42)式をまとめると、(1.30)式は最終的に次のようになる。
𝛼

|𝛼|
𝛽

𝜕𝛹
+ i𝑉𝛹 + (−i𝜵 − 𝑨) 𝛹 − 𝛹 + |𝛹| 𝛹 = 0
𝜕𝑡(1.43)

さらに、(1.43)式の両辺を𝛼(|𝛼|/𝛽) / で除算すると、次のようになる。
𝜕𝛹+ i𝑉𝛹 + (−i𝜵 − 𝑨) 𝛹 − 𝛹 + |𝛹| 𝛹 = 0
𝜕𝑡(1.44)(1.31)式も同様に簡易化する。(1.31)式左辺第一項は、𝑨が一定として次のようになる。

1 |𝛼|𝜈𝛁𝑉 → 𝜈 𝛁 ∗ 𝑉
𝜉 𝑒𝛾
=

=

|𝛼|
𝜈𝛁𝑉
𝜉𝑒 ∗ 𝛾(1.45)

2𝑚∗ |𝛼| |𝛼|∙ ∗ 𝜈𝛁𝑉
ℏ𝑒 𝛾(1.31)式左辺第2項は次のようになる。
11 1
1√2𝜇0 𝐻𝑐
𝛁×𝛁×𝑨→ ∙ 𝛁× 𝛁×
𝑨
𝜇
𝜇 𝜉
𝜉
𝜆
=

√2𝐻
𝛁×𝛁×𝑨
𝜉 𝜆

2𝑚∗ |𝛼| 𝑚∗ |𝛼|=∙ ∗∙ √2𝐻 𝛁 × 𝛁 × 𝑨
ℏ
𝑒 𝜇 𝐻
=

(1.46)(2𝑚∗ |𝛼|)𝛁×𝛁×𝑨ℏ 𝑒 ∗𝜇(1.31)式左辺第3項は次のようになる。
iℏe∗ ∗(𝛹 𝛁𝛹 − 𝛹𝛁𝛹 ∗ )2𝑚∗

11(1.47)

→

iℏe∗
2𝑚∗

=

iℏe∗ |𝛼| 1 ∗∙∙ (𝛹 𝛁𝛹 − 𝛹𝛁𝛹 ∗ )2𝑚∗ 𝛽 𝜉=iℏe∗ |𝛼| 2𝑚∗ |𝛼| ∗(𝛹 𝛁𝛹 − 𝛹𝛁𝛹 ∗ )
∙
∙
2𝑚∗ 𝛽
ℏ

=

i𝑒 ∗ |𝛼| 2𝑚∗ |𝛼| ∗(𝛹 𝛁𝛹 − 𝛹𝛁𝛹 ∗)
2𝑚∗ 𝛽

|𝛼|
𝛽

|𝛼|
1
𝛹∗ 𝛁
𝜉
𝛽

𝛹−

|𝛼|
𝛽

|𝛼|
1
𝛹 𝛁
𝜉
𝛽

𝛹∗

(1.31)式左辺第4項は次のようになる。
𝑒∗
𝑒∗
|𝛹|
𝑨
→
𝑚∗
𝑚∗

|𝛼|
𝛽

𝛹

√2𝜇 𝐻
𝑨
𝜆

=

𝑒 ∗ |𝛼| √2𝜇 𝐻
|𝛹| 𝑨
∙
∙
𝑚∗ 𝛽
𝜆

=

𝑒 ∗ |𝛼|𝑚∗ |𝛼|
|𝛹| 𝑨
∙
∙
𝐻
∙
√2𝜇
𝑚∗ 𝛽
𝑒 ∗𝜇 𝐻=𝑒 ∗ |𝛼| 2𝑚∗|𝛼|
|𝛹| 𝑨
𝑚∗ 𝛽

(1.48)(1.45)式から(1.48)式をまとめると、(1.31)式は最終的に次のようになる。
3(2𝑚∗ |𝛼|)22𝑚∗ |𝛼| |𝛼|∙ ∗ 𝜈𝛁𝑉 +
𝛁×𝛁×𝑨
ℏ
𝑒 𝛾
ℏ2 𝑒∗ 𝜇0+(1.49)i𝑒 ∗ |𝛼| 2𝑚∗ |𝛼| ∗𝑒 ∗ |𝛼| 2𝑚∗ |𝛼|∗)(𝛹|𝛹| 𝑨 = 0
𝛁𝛹
−
𝛹𝛁𝛹
−
2𝑚∗ 𝛽
𝑚∗𝛽

さらに、(1.49)式の両辺を|𝛼| 2𝑚∗ |𝛼|で除算すると、次のようになる。
1
2𝑚∗
𝜈𝛁𝑉
+
𝛁×𝛁×𝑨
ℏ𝑒 ∗ 𝛾ℏ2 𝑒∗𝜇0∗𝑒i+ − ∗ |𝛹| 𝑨 − (𝛹 ∗ 𝛁𝛹 − 𝛹𝛁𝛹 ∗ ) = 0𝑚 𝛽2(1.50)

さらに、(1.50)式を簡単化するために以下の関係を用いる。
𝑱=

1
𝛁×𝛁×𝑨
𝜇

(1.51)

𝛁∙𝑱=0(1.52)

12ここで、J は超伝導体に流れる電流であり、(1.52)式は電流の発散を表す式である。(1.31)式の両辺で𝛁との内積を取ると、(1.51)式、(1.52)式から、次のような式が得られる。
𝛁∙

1
𝑒∗
i
|𝛹| 𝑨 − (𝛹 ∗ 𝛁𝛹 − 𝛹𝛁𝛹 ∗)
𝜈𝛁𝑉
−
∗
∗
ℏ𝑒 𝛾
𝑚 𝛽
2

↔−

𝑚∗ 𝛽
ℏ𝑒 ∗ 𝛾=0i𝜈𝛁 𝑉 = (𝛹 ∗ 𝛁 𝛹 − 𝛹𝛁 𝛹 ∗ ) − 𝛁 ∙ (|𝛹| 𝑨)
2(1.53)i↔ 𝜎𝛁 𝑉 = (𝛹 ∗ 𝛁 𝛹 − 𝛹𝛁 𝛹 ∗ ) − 𝛁 ∙ (|𝛹| 𝑨)
2
ここで、
−

𝑚∗ 𝛽
ℏ𝑒 ∗ 𝛾𝜈→𝜎(1.54)とする。(1.53)式はスカラーポテンシャル𝑉を記述する。本研究の数値解析では、(1.44)式および(1.53)式を解くことによって超伝導体の解析を行う。

1.4 アフィン数値積分法(Affine Integrator(AFI))
アフィン変換数値積分法(AFI)[1][2][6]とは、近年提案された陽的数値積分法である。
拡散方程式のような放物型偏微分方程式や、ゲージ場存在下の Schrödinger 方程式やTime-Dependent Ginzburg-Landau(TDGL) 方程式の数値積分のために考案された。
数値的安定性に関して、AFI はラプラシアン項に関して無条件安定であり、方程式を構成する項全体に関して高い数値的安定性を有する。エネルギー保存に関しては、AFI は線
形保存系(Schrödinger 方程式や TDGL 方程式)において全エネルギーが厳密に保存される性質を持つ。
AFI は空間に関して離散化するための格子をチェッカーボード状に分解することから導かれるアフィン変換対で構成される。数学的に単純なスキームであり、数値的安定性やエ
ネルギー保存に関する理論解析が容易である。また必要とされる記憶領域は最小であり、付加的なワーキングメモリは必要とされないことも AFI の特徴である。

1.4.1

AFI のイメージ

図 1-1 の左側に空間に関する離散化のためのチェッカーボード格子を示す。基礎方程式
はこのチェッカーボード格子によって、まずは空間に関して離散化される。右側にチェッ
カーボード格子に基づいて構成された高次元常微分方程式の概形を示す。
黒丸格子点の集合を Q 集合、白丸格子点の集合を P 集合と呼ぶことにする。Q 集合お
よび P 集合上で標本化された従属変数の時間発展は図 1-1 の右側に示された方程式のブロ13

ック行列が支配する。
図 1-2 に常微分方程式を時間に関して離散化するようすを示す。このように、AFI スキ
ームは Q 集合上で定義された従属変数の時間発展と P 集合上で定義された従属変数の時間発展を交互に繰り返していくものである。

図 1-1

空間に関する離散化のイメージ

図 1-2

時間に関する離散化のイメージ

1.4.2 1 次元空間において AFI 法の時間発展方程式の導出
1 次元空間での AFI の数値積分スキームを導出する。
まずは空間に関する離散化を行う。本研究で用いる基礎方程式は TDGL 方程式である(1.55)式とする。
𝛾𝜕𝜓= (𝛁 − i𝑨) 𝜓 − 𝛼𝜓 − 𝛽|𝜓| 𝜓𝜕𝑡(1.55)

1 次元空間では(1.56)式が(1.55)式で適用されており、ベクトルポテンシャルは固定されていると考える。
𝛁=

𝜕
,
𝜕𝑥

(1.56)

𝑨=𝐴

ベクトルポテンシャルはリンク変数𝑤 , によって次のように実装される。
𝑤 = exp i𝜃,𝜃 = −ℎ𝐴
14

𝑥 +𝑥
,
2

𝑖, 𝑗 = 1,2,3,4.

(1.57)

ここで𝑥 , 𝑥 はそれぞれ格子点𝑖, 𝑗の𝑥座標で、ℎは格子点間隔である。格子点𝑖と格子点𝑗が接続していない場合は𝑤 , = 0と約束する。簡単のため、周期的境界条件を採用する。空間に関する離散化のために導入した格子およびリンク変数の状況を図 1-3 に示す。

1 次元空間での周期的境界条件下の

図 1-3格子点とリンク変数のインデックスの定義(1.55)式を空間に関して離散化する。𝜓 (𝑡) = 𝜓(𝑥 , 𝑡), 𝛼 = 𝛼(𝑥 ), 𝛽 = 𝛽(𝑥 )とおくと(1.58)式で表される連立微分方程式が得られる。
𝑑𝜓
1
⎧𝛾
= (𝑤
ℎ
⎪ 𝑑𝑡
𝑑𝜓
1
⎪𝛾
= (𝑤
𝑑𝑡
ℎ
𝑑𝜓
1
⎨𝛾
= (𝑤
⎪ 𝑑𝑡
ℎ
1
⎪ 𝑑𝜓
⎩ 𝛾 𝑑𝑡 = ℎ ( 𝑤𝜓 + 𝑤 𝜓 − 2𝜓 ) − 𝛼 𝜓 − 𝛽 |𝜓 | 𝜓
𝜓 + 𝑤 𝜓 − 2𝜓 ) − 𝛼 𝜓 − 𝛽 |𝜓 | 𝜓(1.58)

𝜓 + 𝑤 𝜓 − 2𝜓 ) − 𝛼 𝜓 − 𝛽 |𝜓 | 𝜓
𝜓 + 𝑤 𝜓 − 2𝜓 ) − 𝛼 𝜓 − 𝛽 |𝜓 | 𝜓

ここで、𝑤 は𝑤 の複素共役を表す。
奇数の格子点の組(図 1-3 の黒丸に対応)と偶数の格子点の組(図 1-3 の白丸に対応)の順に並べているが、図 1-1 の右側のイメージに従って表記している。この連立微分方程式は以下の(1.59)式で表せる。
𝜓
𝑑 𝜓
𝛾
𝑑𝑡 𝜓
𝜓
−2 − 𝑈 ℎ
1⎛
0
= ⎜
ℎ
𝑤
𝑤
⎝

0
−2 − 𝑈 ℎ𝑤
𝑤𝑤𝑤−2 − 𝑈 ℎ
0

𝑤
⎞
𝑤
⎟
0
−2 − 𝑈 ℎ ⎠

𝜓
𝜓
𝜓
𝜓

(1.59)

ここで、𝑈 =𝛼 +𝛽 𝜓, 𝑖 = 1,2,3,4.

とおく。𝜓 は𝜓 の推定値である。理由は後述する。

15(1.60)

また、
𝜓
𝜓

𝒒=

𝜓
𝜓

,𝒑 =

(1.61)(1.62)

𝜎 = 2 + 𝑈 ℎ , 𝑖 = 1,2,3,4.
とおく。また、
𝐷 −𝜎

=

−𝜎
0

0
−𝜎
𝑤
= 𝑤

, 𝐷(−𝜎 ) =

−𝜎
0

0
−𝜎

,𝑊
(1.63)𝑤𝑤とすると、(1.59)式は(1.64)式で表せる。
𝛾d 𝒒
1 𝐷 −𝜎
=
𝒑
dt
ℎ
𝑊

𝒒
𝒑

𝑊
𝐷(−𝜎 )(1.64)

時間推進演算子 𝐴を用いると(1.64)式は(1.65)式で表される。また、 𝐴は(1.66)で定義さ
れる。
𝛾

𝐴=

1
ℎ

d 𝒒
𝒒
= 𝐴 𝒑
𝒑
dt

−𝜎 𝜓 +∈{ , , , }(1.65)𝑤 𝜓{ , , , }𝜕𝜕𝜓

(1.66)

𝑄 = {1,3}, 𝑃 = {2,4}として時間推進演算子𝐴を次のように分解する。
(1.67)

𝐴=𝐴 +𝐴
ただし、
𝐴 =

𝐴 =

1
ℎ

1
ℎ

−𝜎 𝜓 +
∈

𝑤 𝜓

𝜕
𝜕𝜓

(1.68)𝑤 𝜓𝜕𝜕𝜓

(1.69)∈−𝜎 𝜓 +∈∈である。また、𝑤 = 𝑤 である。
次に、(1.65)式を時間に関して離散化することを考える。
時間刻み幅を𝜏、(𝒒, 𝒑) = (𝒒(𝑡), 𝒑(𝑡))、(𝒒 , 𝒑 ) = (𝒒(𝑡 + 𝜏), 𝒑(𝑡 + 𝜏))として時間発展方程式の(1.70)式が得られる。
𝒒
𝒑

= exp

𝜏
𝐴
𝛾

𝜏
𝒒
𝒑 = exp 𝛾 𝐴 + 𝐴

𝒒𝒑(1.70)

指数関数演算子は Lie-Trotter-Suzuki(LTS)分解によって(1.71)式のように近似できる。
LTS 分解とは、exp 𝑥(𝐴 + 𝐵)の形のものを(1.72)式のように積に分解する方法の 1 つであ
16

る。LTS 分解でどのように分解できるかは後述する。また、𝑂は誤差オーダー部分であ
る。
exp

𝜏
𝐴 +𝐴
𝛾

⎧
⎪

exp

𝜏
𝐴
𝛾

exp

𝜏
𝐴
𝛾

+ 𝑂(𝜏 )

𝜏
⎨exp
𝐴
⎪
2𝛾
⎩

exp

𝜏
𝐴
𝛾

exp

𝜏
𝐴
2𝛾

=

…

𝑥
𝑥
exp 𝑥(𝐴 + 𝐵) = exp 𝐴 exp 𝐵𝑛𝑛
従って、指数関数演算子exp

𝐴

とexp

(1.71)+ 𝑂(𝜏 )

+𝑂𝜏𝑛(1.72)𝐴 の行列表現が得ると、数値積分スキーム(AFI)を構成できる。
まず、ここで演算子𝐴 を𝜓 (𝑖 ∈ 𝑄)に 1 回作用させる。
𝐴 𝜓 =

1
ℎ

−𝜎 𝜓 +

𝑤 𝜓

,

𝑖∈𝑄

(1.73)∈

次に、演算子𝐴 を𝜓 (𝑖 ∈ 𝑄)に 2 回作用させる。
𝐴 𝜓 =

−𝜎
ℎ

𝜓 −

1
𝜎

𝑤 𝜓

,

𝑖∈𝑄

(1.74)∈

従って、演算子𝐴 を𝜓 (𝑖 ∈ 𝑄)に𝑚回(𝑚は正の整数)作用させると
𝐴 𝜓 =

−𝜎
ℎ

𝜓 −

1
𝜎

𝑤 𝜓
∈

となる。

17

,

𝑖 ∈ 𝑄, 𝑚 ≥ 1(1.75)

ゆえに、マクローリン展開をすると
𝜏
exp
𝐴
𝛾

∞

𝜓 =
𝑚=0

𝒎
1 𝜏
𝐴𝑄 𝜓𝑖 = 𝜓𝑖 +𝑚! 𝛾

1 𝜏 𝒎 −𝜎
𝑚! 𝛾
ℎ

=𝜓 +

=𝜓 + 𝜓 −=𝜓 + 𝜓 −

1
𝜎

𝜓 −

𝑚=1

1
𝜎

1 𝜏 𝒎 𝑚
𝐴𝑄 𝜓𝑖
𝑚! 𝛾
𝑤 𝜓

∈

𝑤 𝜓

1 −𝜎 𝜏
𝑚! ℎ 𝛾𝑤 𝜓

1 −𝜎 𝜏
𝑚! ℎ 𝛾

∈

1
𝜎

∞

∈

(1.76)− 𝜓 −=𝜓

1 −𝜎 𝜏
𝑚! ℎ 𝛾

+

1
𝜎

1
𝜎

𝑤 𝜓
∈

𝑤 𝜓

1 −𝜎 𝜏
𝑚! ℎ 𝛾1−∈𝜎𝜏
1 − exp −𝜎𝜏ℎ 𝛾= exp −𝜓 +ℎ 𝛾
𝜎=𝑎 𝜓 +𝑏𝑤 𝜓∈𝑤 𝜓 ,𝑖 ∈ 𝑄∈が得られる。ただし、
𝑎 = exp

−𝜎 𝜏
1−𝑎
,𝑏 =
ℎ 𝛾
𝜎

(1.77)

とおいた。(1.76)式は 1 次元空間における時間発展方程式である。また、𝐴 𝜓 = 0,𝑖∈𝑃

(1.78)𝐴 𝜓 = 0,𝑖∈𝑃(1.79)であるので、であることがわかる。演算子𝐴 でも同様な計算を行うことで次のようになる。
exp

𝜏
𝐴
𝛾

𝒒
𝐷 𝑎
𝒑 =
0

exp

𝜏
𝐴
𝛾

𝒒
0
𝒑 = 𝐷(𝑏 )𝑊

ここで、𝐼は単位行列であり、
18

𝐷 𝑏 𝑊
𝐼

𝒒
𝒑

(1.80)

𝐼
𝐷(𝑎 )

𝒒
𝒑

(1.81)

𝐷 𝑎

𝑎
0

=

0
,𝐷 𝑏
𝑎

𝑎
𝐷(𝑎 ) =
0

𝑏
0

=

0
𝑏
, 𝐷(𝑏 ) =
𝑎
0

0
,
𝑏

(1.82)

0𝑏

とおいた。
(1.80)式と(1.81)式はそれぞれ指数関数演算子exp

𝐴

とexp

𝐴 の行列表現を表しており、数値積分スキーム(AFI)が得られたことになる。
次に、(1.60)式の非線形項𝛽 𝜓で、𝜓 ではなく推定値 𝜓 が用いられる理由を説明する。推定値は(1.83)式から(1.86)式で計算される。
1
𝜓 = (𝑤 𝜓 + 𝑤 𝜓 )2(1.83)1
𝜓 = (𝑤 𝜓 + 𝑤 𝜓 )2(1.84)1
𝜓 = (𝑤 𝜓 + 𝑤 𝜓 )2(1.85)1
𝜓 = (𝑤 𝜓 + 𝑤 𝜓 )2(1.86)

注意する点として、𝜓 (𝑖 ∈ 𝑄)の値は𝜓 (𝑖 ∈ 𝑃)で推定され、𝜓 (𝑖 ∈ 𝑃)の値は𝜓 (𝑖 ∈ 𝑄)で推定されることである。このことにより、指数関数演算子の行列表現構成が容易になる。
例を挙げると、指数関数演算子exp

𝐴の行列表現を得るために、(1.75)式において𝐴をベクトル𝒒の要素に何回も作用させた。その際、𝜎 , 𝑖 ∈ 𝑄はベクトル𝒒の要素に依存せず、ベクトル𝒑の要素に依存することが重要である。つまり、時間推進演算子𝐴 にとって𝜎 , 𝑖 ∈ 𝑄は定数となる。このことが行列表現の導出を可能としている。

1.4.3 2 次元空間において AFI 法の時間発展方程式の導出
2 次元空間での AFI 法の数値積分スキームを導出する。
基礎方程式は 1 次元空間で用いたものと同様に(1.55)式を用いる。ただし、
𝛁=

𝜕 𝜕
,
,𝑨 = 𝐴 ,𝐴𝜕𝑥 𝜕𝑦(1.87)である。
ここで、格子点やリンク変数を指定するためにダブルインデックスを用いる。そうする
ことでダブルインデックスを格子点やリンク変数を空間的な位置に対応させることができ19る。図 1-4 にダブルインデックスの定義を示す。また、リンク変数は𝑥成分の𝑤 , と、𝑦成分の𝑤 , に分けて定義される。境界条件は𝑥方向、𝑦方向ともに周期的境界条件とする。

2 次元空間での周期的境界条件下の

図 1-4格子点とリンク変数のダブルインデックスでの定義

まず空間に関して離散化すると(1.88)式が得られる。
γ

𝑑𝜓 ,
𝑑𝑡

=

1
𝑦
𝑦
𝑤𝑥𝑖,𝑗 𝜓𝑖+1,𝑗 + 𝑤𝑥𝑖−1,𝑗 𝜓𝑖−1,𝑗 + 𝑤𝑖,𝑗 𝜓𝑖,𝑗+1 + 𝑤𝑖,𝑗−1 𝜓𝑖,𝑗−1 − 4𝜓 ,ℎ−𝛼 , 𝜓 , − 𝛽 , 𝜓 ,(1.88)

𝜓,

ただし、
𝑈, = 𝛼 +𝛽, 𝜓 ,
𝜓, =

1 𝑥
𝑦
𝑦
𝑤 𝜓
+ 𝑤𝑥𝑖−1,𝑗 𝜓𝑖−1,𝑗 + 𝑤𝑖,𝑗 𝜓𝑖,𝑗+1 + 𝑤𝑖,𝑗−1 𝜓𝑖,𝑗−1
4 𝑖,𝑗 𝑖+1,𝑗(1.89)(1.90)

とおく。また、𝜎, =4+𝑈, ℎ(1.91)

とおく。
ここまで格子点とリンク変数をダブルインデックスで定義したが、図 1-4 において格子点を(1,1)、(1,2)、…、(4,3)、(4,4)と定義している点をシングルインデックス 1、
２、…、15、16 と定義したものを考えると、1 次元空間で導出した数値積分スキームの20

(1.80)式と(1.81)式、時間発展方程式の(1.76)式がそのまま用いることができる。数値積分
スキームは図 1-4 の黒丸に対して(1.80)式、白丸に対して(1.81)式が対応する。
具体例として格子点(1,1)での時間発展方程式を(1.92)式に示す。なお、従属変数𝜓 , が時間発展した後のものを𝜓′ , とする。(1.76)式を適用するが、使用されている格子点とリンク変数に関してはダブルインデックスで定義した表記を使用している。ほかの格子点についても同様に(1.76)式を用いて導出可能である。
𝜓

𝑦

,

𝑦

= 𝑎 , 𝜓 , + 𝑏 , 𝑤𝑥1,1 𝜓2,1 + 𝑤𝑥4,1 𝜓4,1 + 𝑤1,1 𝜓1,2 + 𝑤1,4 𝜓1,4
𝑎 , = exp −𝜎 , 𝜏⁄ℎ(1.93)𝑏 , = 1−𝑎 ,(1.94)𝜎,𝜎 , =4+ 𝛼 , +𝛽 , 𝜓 ,𝜓, =(1.92)(1.95)ℎ

1(𝑤 , 𝜓 , + 𝑤 , 𝜓 , + 𝑤 , 𝜓 , + 𝑤 , 𝜓 , )4(1.96)

これを一般化すると、2 次元空間における AFI 法での時間発展方程式の(1.97)式を得る
ことができる。
𝜓 , =𝑎, 𝜓,+𝑏 , 𝑤 , 𝜓

,

+𝑤

, 𝜓

,

+𝑤, 𝜓 ,

1
𝑤, 𝜓
4

,

𝜓,

(1.97)

𝑎 , = exp −𝜎 , 𝜏⁄ℎ(1.98)𝑏, = 1−𝑎,(1.99)𝜎,𝜎, = 4+ 𝛼 , +𝛽, 𝜓 ,
𝜓, =

+𝑤,

+𝑤

, 𝜓

,

+𝑤, 𝜓,(1.100)ℎ+𝑤,

𝜓,(1.101)

1.5 Lie-Trotter-Suzuki(LTS)分解
ここでは、指数関数演算子
exp 𝜏(𝐴 + 𝐵)

(1.102)の Lie-Trotter-Suzuki(LTS)分解[3][4]について考える。LTS 分解は、(1.102)式を積に分解する方法の一つである。一般的に、演算子𝐴と𝐵は非可換である。
一般的な LTS 分解は次のようになる。
exp 𝜏(𝐴 + 𝐵) = exp(𝜏𝐴) exp(𝜏𝐵) + 𝑂(𝜏 )(1.103)

AFI における 2 次のオーダーの LTS 分解として次の式を用いる。
exp 𝜏(𝐴 + 𝐵) = exp (1/2)𝜏𝐴 exp(𝜏𝐵) exp (1/2)𝜏𝐴 + 𝑂(𝜏 )(1.104)しかし、近年最適な LTS 分解が提案された。(1.105)式は 2 次のオーダーの LTS 分解として用いられるが、さらに最適化された 2 次のオーダーの LTS 分解は次のようになる。

21

exp 𝜏(𝐴 + 𝐵)
= exp(𝑎 𝜏𝐴) exp(𝑏 𝜏𝐵) exp(𝑎 𝜏𝐴) exp(𝑏 𝜏𝐵) exp(𝑎 𝜏𝐴)(1.105)+ 𝑂(𝜏 )𝑎 =1
3 − √3 ,6

𝑎 =1−𝑎 ,𝑏 =1
2

また、最適化された 4 次のオーダーの LTS 分解は次のようになる。
exp 𝜏(𝐴 + 𝐵)
= exp(𝑎 𝜏𝐴) exp(𝑏 𝜏𝐵) exp(𝑎 𝜏𝐴) exp(𝑏 𝜏𝐵) exp(𝑎 𝜏𝐴) exp(𝑏 𝜏𝐵)(1.106)× exp(𝑎 𝜏𝐴) exp(𝑏 𝜏𝐵) exp(𝑎 𝜏𝐴) exp(𝑏 𝜏𝐵) exp(𝑎 𝜏𝐴) + 𝑂(𝜏 )𝑎 = 0.095848502741203681182,
𝑎 =0.078111158921637922695,
𝑎 =

1
2

(𝑎 + 𝑎 ),𝑏 = 0.42652466131587616168,
𝑏 =0.12039526945509726545,
𝑏 =1

2(𝑏 + 𝑏 ),さらに最適化された 6 次のオーダーの LTS 分解は次のようになる。
exp 𝜏(𝐴 + 𝐵) =exp

1
1
𝑤 𝜏𝐴 exp(𝑤 𝜏𝐵) exp (𝑤 + 𝑤 )𝜏𝐴
2
2

× exp(𝑤 𝜏𝐵) exp
× exp

1(𝑤 + 𝑤 )𝜏𝐴 exp(𝑤 𝜏𝐵)
2

1
1(𝑤 + 𝑤 )𝜏𝐴 exp(𝑤 𝜏𝐵) exp (𝑤 + 𝑤 )𝜏𝐴
2
2

1
× exp(𝑤 𝜏𝐵) exp (𝑤 + 𝑤 )𝜏𝐴 exp(𝑤 𝜏𝐵)
2
× exp

1
1
(𝑤 + 𝑤 )𝜏𝐴 exp(𝑤 𝜏𝐵) exp (𝑤 + 𝑤 )𝜏𝐴
2
2

× exp(𝑤 𝜏𝐵) exp

11(𝑤 + 𝑤 )𝜏𝐴 exp(𝑤 𝜏𝐵) exp 𝑤 𝜏𝐴
2
2

+𝑂(𝜏 )
𝑤 = 0.18793069262651671457
𝑤 = 0.5553
22(1.107)

𝑤 = 0.12837035888423653774
𝑤 = −0.84315275357471264676
𝑤 = 1 − 2(𝑤 + 𝑤 + 𝑤 + 𝑤 )1.6 研究の目的
常微分方程式を数値解析する際に Euler 法や Runge-Kutta 法がよく用いられる。Euler 法や Runge-Kutta 法は常微分方程式の数値解法で、数学的に簡単でありプログラムを記述する上でも簡単である。
対して AFI は、Euler 法や Runge-Kutta 法と同じく数学的に簡単であり、プログラムを記述する上でも簡単である。さらに、方程式を構成する項全体に関して高い数値的安定性を有している。したがって、拡散方程式や量子拡散方程式、TDGL 方程式の数値解析を行う上で AFI 法を用いると高い安定性のもと数値解析を行うことができる。
本研究では、拡散方程式、量子拡散方程式、TDGL 方程式をシミュレーションプログラ
ムを作成する。その時に用いる数値解法を Euler 法、Runge-Kutta 法、AFI と変え、時間刻
み幅を変更しながら誤差評価や TDGL 方程式の AFI を用いた時の安定性を調査することを目的とする。

23

第 2 章 実装方法と計算条件
2.1 実装方法
実装するにあたって、統合開発環境でありプログラミング言語でもある Processing で実装を行う。Processing は、アーティストよるコンテンツ制作作業のために、詳細な設定を行
う関数を排除しており、非プログラマーでも容易にプログラムを実装できる。
シミュレーション領域は 1 次元では 1 次元直線であり、直線上に超伝導領域が存在するとしてプログラムを実装する。また、シミュレーション領域は 2 次元では 2 次元平面であり、平面上に超伝導領域が存在するとしてプログラムを実装する。
時間刻み幅を変更して Euler 法や Runge-Kutta 法、AFI を用いて得られる数値解と解析解から得られる誤差を一定時間計測する。その計測する間で得られる最大の誤差を、使用した時間刻み幅の数値解を得るのに使用した時間刻み幅の用いた数値解法における誤差として評価する。
またはある時間刻み幅とその半分の時間刻み幅を用いてから数値解を計測し、それを用いて誤差推定を行うことで使用した時間刻み幅の用いた数値解法における誤差として評価する。1 次元 TDGL 方程式の安定性は、非線形のため誤差推定から算出する。

2.2 1 次元拡散方程式の誤差評価
ここでは数値計算の精度の評価[6]のために 1 次元拡散方程式を実装した。定義域は0 ≤𝑥 ≤ 𝐿であるとする。従属変数𝜓(𝑥, 𝑡)が従う基礎方程式は次のようになる。また、周期的境界条件を採用する。
𝜕𝜓 𝜕 𝜓
=
𝜕𝑡
𝜕𝑥

(2.1)

この式を数値的に解いて得られた𝜓の値と解析解から得られる𝜓の値の差を用いて計算される誤差を使用した。ただし、格子点間隔をℎ、格子点の数を𝑁、定義域の⾧さを𝐿とする。
定義域の⾧さ𝐿と関係は以下のようになる。
(2.2)

𝐿 = 𝑁ℎここで、周期的境界条件として、以下のようにする。
𝜓(0, 𝑡) = 𝜓(𝐿, 𝑡)(2.3)

𝜓(𝐿 + ℎ, 𝑡) = 𝜓(1, 𝑡)(2.4)

初期条件は次のようにする。
𝜓 (0) = 𝜓(𝑖ℎ, 0) = sin(𝑘𝑖ℎ) ,𝑘=

まず、解析解の式を表すと以下のようになる。
24

2𝜋𝑚
,
𝐿

𝑚 = 1,2,3, ⋯(2.5)

𝜓(𝑥, 𝑡) = exp(−𝑘 𝑡) sin(𝑘𝑥)(2.6)

次に、数値解について説明する。
Euler 法においては、従属変数𝜓 は次の式に従って、時間𝜏後の値𝜓 に更新される。
𝜓 =𝜓 +𝜏(𝜓ℎ+𝜓− 2𝜓 )(2.7)

ここで時間𝑡を以下の式が成り立つようにと置くと、
𝑡 = 𝑛𝜏,𝑛 = 0,1,2, ⋯(2.8)

その時間における誤差𝐸𝑟𝑟(𝑡)、その時間刻み幅𝜏における誤差𝐸𝑅𝑅は以下のようになる。ここで𝑇は十分な大きさの定数である。
𝐸𝑟𝑟(𝑡) =1𝑁𝜓 (𝑡) − 𝜓(𝑥, 𝑡)

𝐸𝑅𝑅 = max 𝐸𝑟𝑟(𝑡)(2.9)(2.10)

2 次の Runge-Kutta 法(RK2)、4 次の Runge-Kutta 法(RK4)、2 次のオーダーの LTS 分解を使用した AFI(AFI-2)、最適化された 2 次のオーダーの LTS 分解を使用した AFI(AFI-2opt)、最適化された 4 次のオーダーの LTS 分解を使用した AFI(AFI-4-opt)、最適化された
6 次のオーダーの LTS 分解を使用した AFI(AFI-6-opt)の場合でもその時間刻み幅𝜏における誤差𝐸𝑅𝑅を算出する。
図 2-1 に(2.1)式に基づく 1 次元拡散方程式の数値シミュレーションの実行例を示す。白棒が解析解、黒丸が数値解を示す。初期条件として𝜓 (0) = sin𝑖ℎ 、境界条件として周期的境界条件を用いた。ただし、𝑚 = 1、ℎ = 0.1、𝑁 = 100、𝐿 = 10、𝜏 = 0.005である。

図2-1

1 次元拡散方程式のシミュレーション
25

2.3 2 次元拡散方程式の誤差評価
ここでは数値計算の精度の評価のために 2 次元拡散方程式を実装した。定義域は正方形の𝑥 ∈ [0, 𝐿] × 𝑦 ∈ [0, 𝐿]であるとする。従属変数𝜓(𝑥, 𝑦, 𝑡)が従う基礎方程式は次のようになる。
また、周期的境界条件を採用する。
𝜕𝜓
𝜕
𝜕
=
+
𝜕𝑡
𝜕𝑥
𝜕𝑦

(2.11)

𝜓

この式を数値的に解いて得られた𝜓の値と解析解から得られる𝜓の値の差を用いて計算される誤差を使用した。ただし、格子点間隔をℎ、格子点の数を𝑁、定義域の⾧さを𝐿とする。
定義域の⾧さ𝐿と関係は以下のようになる。
(2.12)

𝐿 = 𝑁ℎここで、周期的境界条件として、以下のようにする。
𝑥方向(左端と右端)は次のようになる。

(2.13)𝑦 ∈ [0, 𝐿]

𝜓(0, 𝑦, 𝑡) = 𝜓(𝐿, 𝑦, 𝑡),(2.14)𝑦 ∈ [0, 𝐿]

𝜓(𝐿 + ℎ, 𝑦, 𝑡) = 𝜓(1, 𝑦, 𝑡),𝑦方向(上端と下端)は次のようになる。

(2.15)𝑥 ∈ [0, 𝐿]

𝜓(𝑥, 0, 𝑡) = 𝜓(𝑥, 𝐿, 𝑡),(2.16)𝑥 ∈ [0, 𝐿]

𝜓(𝑥, 𝐿 + ℎ, 𝑡) = 𝜓(𝑥, 1, 𝑡),
初期条件は次のようにする。

𝜓 , (0) = 𝜓(𝑖ℎ, 𝑗ℎ, 0) = sin(𝑘 𝑖ℎ) sin 𝑘 𝑗ℎ ,
𝑘 =

2𝜋𝑚
2𝜋𝑛
,𝑘 =
,
𝐿
𝐿

𝑚, 𝑛 = 1,2,3, ⋯(2.17)

まず、解析解の式を表すと以下のようになる。
𝜓(𝑥, 𝑦, 𝑡) = exp − 𝑘 + 𝑘 𝑘 𝑡 sin(𝑘 𝑥) sin 𝑘 𝑦(2.18)

次に、数値解について説明する。
Euler 法においては、従属変数𝜓 , は次の式に従って、時間𝜏後の値𝜓 , に更新される。
𝜓, =𝜓, +

𝜏
𝜓
ℎ

,

+𝜓

,

+𝜓 ,

+𝜓,

− 4𝜓 ,(2.19)

ここで時間𝑡を(2.8)式が成り立つようにと置くと、その時間における誤差𝐸𝑟𝑟(𝑡)、その時
間刻み幅𝜏における誤差𝐸𝑅𝑅は以下のようになる。ここで𝑇は十分な大きさの定数である。
𝐸𝑟𝑟(𝑡) =1𝑁𝜓 , (𝑡) − 𝜓(𝑥, 𝑦, 𝑡)

𝐸𝑅𝑅 = max 𝐸𝑟𝑟(𝑡)(2.20)(2.21)

RK2、RK4、AFI-2、AFI-2-opt、AFI-4-opt の場合でもその時間刻み幅𝜏における誤差
𝐸𝑅𝑅を算出する。
26

図 2-2 に(2.11)式に基づく 2 次元拡散方程式の数値シミュレーションの実行例を示す。左側が解析解、右側が数値解を示す。初期条件として𝜓 , (0) = sin𝑖ℎ sin𝑗ℎ 、境界条件

として周期的境界条件を用いた。ただし、𝑚 = 1、𝑛 = 1、ℎ = 0.1、𝑁 = 100、𝐿 = 10、𝜏 =0.001である。

図2-2

2 次元拡散方程式のシミュレーション

2.4 1 次元量子拡散方程式の誤差評価(ディリクレ境界条件)
ここでは数値計算の精度の評価のために 1 次元量子拡散方程式を実装した。定義域は0 ≤𝑥 ≤ 𝐿であるとする。従属変数𝜓(𝑥, 𝑡)が従う基礎方程式は次のようになる。また、ディリクレ境界条件を採用する。
𝛾

𝜕𝜓
𝜕
=
−𝐴
𝜕𝑡
𝜕𝑥

𝜓

(2.22)

この式を数値的に解いて得られた𝜓の値と解析解から得られる𝜓の値の差を用いて計算される誤差を使用した。ただし、定義域の⾧さを𝐿、ベクトルポテンシャルの𝑥成分𝐴 が一様であるとする。格子点間隔をℎ、格子点の数を𝑁、定義域の⾧さを𝐿とする。定義域の⾧さ𝐿と関係は以下のようになる。
𝐿 = 𝑁ℎ(2.23)

ここで、ディリクレ境界条件として、以下のようにする。
𝜓(0, 𝑡) = 𝜓(𝐿 + ℎ, 𝑡) = 0
初期条件は次のようにする。

27(2.24)

𝜓 (0) = 𝜓(𝑖ℎ, 0) = sin(𝑘𝑖ℎ) ,

𝑘=

𝜋𝑚
,
𝐿+ℎ

𝑚 = 1,2,3, ⋯(2.25)

時定数𝛾を次のようにする。
𝛾(2.26)= 𝜂 + i𝜔

リンク変数wを次のようにする。
(2.27)

𝑤 = exp(−iℎ𝐴 )ここで、𝑤は𝑤の複素共役を表す。
まず、解析解の式を表すと以下のようになる。
𝜓(𝑥, 𝑡) =1 exp[−𝜂(𝑘 − 𝐴 ) 𝑡] exp[i(𝑘𝑥 − 𝜔(𝑘 − 𝐴 ) 𝑡)]
2i − exp[−𝜂(𝑘 + 𝐴 ) 𝑡] exp[i(𝑘𝑥 + 𝜔(𝑘 + 𝐴 ) 𝑡)](2.28)

次に、数値解について説明する。
Euler 法においては、従属変数𝜓 は次の式に従って、時間𝜏後の値𝜓 に更新される。
𝜓 =𝜓 +

𝜏
(𝑤𝜓
ℎ

+ 𝑤𝜓

− 2𝜓 )(2.29)

ここで時間𝑡を(2.8)式が成り立つようにと置くと、その時間における誤差𝐸𝑟𝑟(𝑡)、その時
間刻み幅𝜏における誤差𝐸𝑅𝑅はそれぞれ(2.9)式、(2.10)式から求められる。
RK2、RK4、AFI-2、AFI-2-opt、AFI-4-opt の場合でもその時間刻み幅𝜏における誤差𝐸𝑅𝑅を算出する。

28

図 2-3 に(2.22)式に基づく 1 次元量子拡散方程式の数値シミュレーションの実行例を示
す。白棒が解析解、黒丸が数値解を示す。初期条件として𝜓 (0) = sin𝑖ℎ 、境界条件としてディリクレ境界条件を用いた。ただし、𝑚 = 1、ℎ = 0.1、𝑁 = 100、𝐿 = 10、𝐴 = 1.2、
𝜏 = 0.005である。

図2-3

1 次元量子拡散方程式のシミュレーション(ディリクレ境界条件)

2.5 1 次元量子拡散方程式の誤差評価(周期的境界条件)
ここでは数値計算の精度の評価のために 1 次元量子拡散方程式を実装した。定義域は0 ≤𝑥 ≤ 𝐿であるとする。従属変数𝜓(𝑥, 𝑡)が従う基礎方程式は次のようになる。また、周期的境界条件を採用する。
𝛾

𝜕𝜓
𝜕
=
−𝐴
𝜕𝑡
𝜕𝑥

𝜓

(2.30)

この式を数値的に解いて得られた𝜓の値と解析解から得られる𝜓の値の差を用いて計算される誤差を使用した。ただし、定義域の⾧さを𝐿、ベクトルポテンシャルの𝑥成分𝐴 が一様であるとする。格子点間隔をℎ、格子点の数を𝑁、定義域の⾧さを𝐿とする。定義域の⾧さ𝐿と関係は以下のようになる。
𝐿 = 𝑁ℎここで、周期的境界条件として設定し、以下のようにする。
29(2.31)

𝜓(0, 𝑡) = 𝜓(𝐿, 𝑡)(2.32)

𝜓(𝐿 + ℎ, 𝑡) = 𝜓(1, 𝑡)(2.33)

初期条件は次のようにする。
2𝜋𝑚,
𝐿𝑚 = 1,2,3, ⋯(2.34)

1 exp[−𝜂(𝑘 − 𝐴 ) 𝑡] exp[i(𝑘𝑥 − 𝜔(𝑘 − 𝐴 ) 𝑡)]
2i − exp[−𝜂(𝑘 + 𝐴 ) 𝑡] exp[i(𝑘𝑥 + 𝜔(𝑘 + 𝐴 ) 𝑡)](2.35)

𝜓 (0) = 𝜓(𝑖ℎ, 0) = sin(𝑘𝑖ℎ) ,𝑘=

まず、解析解の式を表すと以下のようになる。
𝜓(𝑥, 𝑡) =次に、数値解について説明する。
Euler 法においては、従属変数𝜓 は(2.29)式に従って、時間𝜏後の値𝜓 に更新される。
ここで時間𝑡を(2.8)式が成り立つようにと置くと、その時間における誤差𝐸𝑟𝑟(𝑡)、その時
間刻み幅𝜏における誤差𝐸𝑅𝑅はそれぞれ(2.9)式、(2.10)式から求められる。
RK2、RK4、AFI-2、AFI-2-opt、AFI-4-opt の場合でもその時間刻み幅𝜏における誤差𝐸𝑅𝑅を算出する。
図 2-4 に(2.30)式に基づく 1 次元量子拡散方程式の数値シミュレーションの実行例を示
す。白棒が解析解、黒丸が数値解を示す。初期条件として𝜓 (0) = sin𝑖ℎ 、境界条件として周期的境界条件を用いた。ただし、𝑚 = 1、ℎ = 0.1、𝑁 = 100、𝐿 = 10、𝐴 = 1.2、𝜏 = 0.005
である。

図2-4

1 次元量子拡散方程式のシミュレーション(周期的境界条件)30

2.6 1 次元 TDGL 方程式の誤差評価
ここでは数値計算の精度の評価のために 1 次元 TDGL 方程式を実装した。定義域は0 ≤𝑥 ≤ 𝐿であるとする。従属変数𝜓(𝑥, 𝑡)が従う基礎方程式は次のようになる。また、周期的境界条件を採用する。
𝛾

𝜕𝜓
𝜕
=
−𝐴
𝜕𝑡
𝜕𝑥

𝜓 − 𝛼𝜓 − 𝛽|𝜓| 𝜓(2.36)

この式を時間刻み幅𝜏と で数値的に解いて得られた𝜓の値との差を用いて計算される誤差を使用した。ただし、定義域の⾧さを𝐿、ベクトルポテンシャルの𝑥成分𝐴 が一様であるとする。格子点間隔をℎ、格子点の数を𝑁、定義域の⾧さを𝐿とする。定義域の⾧さ𝐿と関係は以下のようになる。
(2.37)

𝐿 = 𝑁ℎここで、周期的境界条件として設定し、以下のようにする。
𝜓(0, 𝑡) = 𝜓(𝐿, 𝑡)(2.38)

𝜓(𝐿 + ℎ, 𝑡) = 𝜓(1, 𝑡)(2.39)

初期条件は次のようにする。
𝜓 (0) = 𝜓(𝑖ℎ, 0) = sin(𝑘𝑖ℎ) ,

𝑘=

2𝜋𝑚
,
𝐿

𝑚 = 1,2,3, ⋯(2.40)

時間刻み幅𝜏の数値解の時間𝑡における空間平均は次のようになる。
𝑅 (𝑡) =1𝑁|𝜓 (𝑡)|(2.41)

誤差推定 [5]は次のように定義される。
∆𝑅(𝜏, 𝑡) = 𝑅 (𝑡) − 𝑅 / (𝑡)(2.42)

AFI-2、AFI-2-opt、AFI-4-opt を使用して数値計算を行い、誤差推定∆𝑅を算出する。

31

図 2-5 に(2.36)式に基づく 1 次元 TDGL 方程式の数値シミュレーションの実行例を示す。
水色の棒が時間刻み幅𝜏の数値解、赤丸が時間刻み幅𝜏/2の数値解を示す。初期条件として
𝜓 (0) = sin𝑖ℎ 、境界条件として周期的境界条件を用いた。ただし、𝑚 = 1、ℎ = 0.1、𝑁 =100、𝐿 = 10、𝐴 = 0.0、𝛼 = −20、𝛽 = 20、𝜏 = 0.01である。

図2-5

1 次元 TDGL 方程式のシミュレーション

2.7 1 次元 TDGL 方程式の安定性
安定性の評価のために 1 次元 TDGL 方程式を実装した。定義域は0 ≤ 𝑥 ≤ 𝐿であるとする。従属変数𝜓(𝑥, 𝑡)が従う基礎方程式は次のようになる。また、周期的境界条件を採用する。
𝛾

𝜕𝜓
𝜕
=
−𝐴
𝜕𝑡
𝜕𝑥

𝜓 − 𝛼𝜓 − 𝛽|𝜓| 𝜓(2.43)

時間刻み幅𝜏と展開係数𝛽(𝛽 = −𝛼)を変更したときの数値的安定性から評価を行った。安定か不安定かの判定は、(2.42)式のように誤差推定を行い、時間刻み幅𝜏を大きくしていっても∆𝑅の値が滑らかに変化している時間刻み幅𝜏までを安定、∆𝑅の値が急激に増加した時
間刻み幅𝜏からを不安定とした。

32

第 3 章 実装結果および考察
3.1 1 次元拡散方程式の誤差評価
𝑚 = 1, 𝐿 = 10, 𝑁 = 100, ℎ = 0.1, 𝑇 = 10における周期的境界条件の 1 次元拡散方程式の誤差の時間刻み幅依存性を図 3-1 に示す。横軸は時間刻み幅𝜏、縦軸は誤差𝐸𝑅𝑅である。

図3-1

1 次元拡散方程式の誤差の時間刻み幅依存性

時間刻み幅に関わらず、Euler 法、RK2 法、RK4 法の方が AFI-2、AFI-2-opt に比べて誤差が小さい。しかし、AFI-4-opt や AFI-6-opt は Euler 法、RK2 法、RK4 法と誤差の大きさはあまり変わらなかった。Euler 法が𝜏 = 2.0 × 10 付近、AFI-4-opt が𝜏 = 6.0 × 10 付近で一時的に誤差の大きさが小さくなっているが、これは数値シミュレーションを実行するときに行われる数値計算で生まれる誤差や丸め誤差によって誤差が小さくなっているだけである。
また、時間刻み幅を大きくしていくと、Euler 法、RK2 法は時間刻み幅𝜏 = 5.0 × 10 まで、RK4 法は𝜏 = 6.0 × 10 までしか数値シミュレーションが行えない。これに比べて AFIは AFI-2、AFI-2-opt は計測した最大の時間刻み幅𝜏 = 1.0 × 10 まで数値シミュレーション
が行え、AFI-4-opt では𝜏 = 6.0 × 10 まで安定して数値シミュレーションが行えた。AFI-633

opt も同じくらいの時間刻み幅の大きさまで数値シミュレーションが行えると予想し実行
したが、𝜏 = 1.0 × 10 までしか安定して数値シミュレーションが行えなかった。このように AFI-4-opt や AFI-6-opt が AFI-2、AFI-2-opt と比較して小さな時間刻み幅で数値的安定性を失うのは、AFI の特性として LTS 分解に負の数を使用しているからである。
これから、AFI は Euler 法、RK2 法、RK4 法といった既存の方法よりも時間刻み幅を大きくしても数値的安定性が大きい AFI-2、AFI-2-opt、AFI-4-opt について考えることにする。

3.2 2 次元拡散方程式の誤差評価
𝑚 = 1, 𝑛 = 1, 𝐿 = 10, 𝑁 = 100, ℎ = 0.1, 𝑇 = 10における周期的境界条件の 2 次元拡散方程式の誤差の時間刻み幅依存性を図 3-2 に示す。横軸は時間刻み幅𝜏、縦軸は誤差𝐸𝑅𝑅である。

図3-2

2 次元拡散方程式の誤差の時間刻み幅依存性

1 次元拡散方程式の誤差評価である図 3-1 と比較する。1 次元から 2 次元に拡張したため誤差は全体的に大きくなっていたが、1 次元拡散方程式と 2 次元拡散方程式の誤差の時間刻み幅依存性にあまり変化は見られない。しかし、Euler 法、RK2 法は𝜏 = 2.0 × 10 、RK4 法
は𝜏 = 3.0 × 10 、AFI-4-opt は𝜏 = 3.0 × 10 までしか安定して数値シミュレーションが行ず、1 次元拡散方程式の時に比べて小さな時間刻み幅でしか安定して数値シミュレーションが行えない。これは、1 次元の時に比べて誤差が大きくなったからである。
しかし、図 3-2 のように 2 次元拡散方程式の誤差の時間刻み幅依存性をみても Euler 法、
34

RK2 法、RK4 法と比較して AFI は時間刻み幅が大きくても安定して数値シミュレーションが行える点が優れている。

3.3 1 次元量子拡散方程式の誤差評価(ディリクレ境界条件)
𝑚 = 1, 𝐿 = 10, 𝑁 = 100, ℎ = 0.1, 𝛾 = 1, 𝑇 = 10におけるディリクレ境界条件の 1 次元量子拡散方程式の誤差の時間刻み幅依存性を図 3.3 に示す。横軸は時間刻み幅𝜏、縦軸は誤差𝐸𝑅𝑅である。1 次元量子拡散方程式の𝐴 は丸で𝐴 = 0、三角で𝐴 = 1.2、四角で𝐴 = 2.4で示し
た。

図3-3

ディリクレ境界条件の 1 次元量子拡散方程式の誤差の時間刻み幅依存性

1 次元量子拡散方程式の𝐴 = 0の誤差の時間刻み幅依存性、つまり丸のグラフは 1 次元拡散方程式の誤差の時間刻み幅依存性は図 3-1 と一致するはずである。図 3-1 のグラフと図3-3 の丸のグラフを比較すると、境界条件が違うので多少変化があるがおおむね一致することが分かる。
ベクトルポテンシャルがあるとき、時間刻み幅に関わらず誤差が大きく、また Euler 法やAFI といった方法に関係なく誤差に違いが見られなかった。これはベクトルポテンシャルがあるとき、境界条件がディリクレ境界条件のため境界付近で解析解と数値解の振る舞いが違く、誤差が大きくなっている。
35

しかし、ベクトルポテンシャルがあるなしに関わらず、Euler 法、RK2 法は時間刻み幅𝜏 =5.0 × 10 まで、RK4 法は𝜏 = 6.0 × 10 まで、AFI-4-opt は𝜏 = 3.0 × 10 までしか数値シ
ミュレーションが行えない。これよりベクトルポテンシャルに関わらず、AFI は Euler 法やRunge-Kutta 法よりも大きい時間刻み幅でも安定してシミュレーションを行うことができる。

3.4 1 次元量子拡散方程式の誤差評価(周期的境界条件)
𝑚 = 1, 𝐿 = 10, 𝑁 = 100, ℎ = 0.1, 𝛾 = 1, 𝑇 = 10における周期的境界条件の 1 次元量子拡散方程式の誤差の時間刻み幅依存性を図 3-4 に示す。横軸は時間刻み幅𝜏、縦軸は誤差𝐸𝑅𝑅である。1 次元量子拡散方程式の𝐴 は丸で𝐴 = 0、三角で𝐴 = 1.2、四角で𝐴 = 2.4で示した。

図3-4

周期的境界条件の 1 次元量子拡散方程式の誤差の時間刻み幅依存性

1 次元量子拡散方程式の𝐴 = 0の誤差の時間刻み幅依存性、つまり丸のグラフは 1 次元拡散方程式の誤差の時間刻み幅依存性は図 3-1 と一致するはずである。図 3-1 のグラフと図3-3 の丸のグラフを比較すると、おおむね一致することが分かる。
ベクトルポテンシャルがあるとき、ベクトルポテンシャルが大きい方が誤差が大きくなるが、安定して数値シミュレーションが行える時間刻み幅はベクトルポテンシャルがない
ときと変わらなかった。
36

図 3-3 と図 3-4 を比べると、ベクトルポテンシャルがないとき、ディリクレ境界条件と周期的境界条件の誤差に違いは見られない。しかし、ベクトルポテンシャルがあるとき、周期的境界条件の方がディリクレ境界条件よりも誤差が小さい。このことは、解析解と数値解の振る舞いの影響である。ベクトルポテンシャルがないときはディリクレ境界条件で得られる数値解が解析解と一致する振る舞いをしているが、図 2-3、図 2-4 に見られるようにベク
トルポテンシャルがあるときは周期的境界条件で得られる数値解が解析解と一致する振る舞いをしているからである。
これからベクトルポテンシャルに関わらず、AFI は Euler 法や Runge-Kutta 法に比べても誤差の大きさは変わらず、また大きな時間刻み幅でも安定して数値計算が行えることが分かる。

3.5 1 次元 TDGL 方程式の誤差評価
𝑚 = 1, 𝐿 = 10, 𝑁 = 100, ℎ = 0.1, 𝛾 = 1, 𝐴 = 0, 𝛼 = −20, 𝛽 = 20における時間刻み幅𝜏を変更したときの周期的境界条件の 1 次元 TDGL 方程式の誤差の時間依存性を図 3-5 に示す。このとき AFI-2 を使用して数値計算を行った。横軸は時刻𝑡、縦軸は誤差推定∆𝑅である。

図3-5

1 次元 TDGL 方程式の誤差の時間依存性(AFI-2)

時刻𝑡を変化させたとき、誤差推定∆𝑅が変化している図 3-5 の左側の 5 点を抜き出して考37

える。抜き出した 5 点を横軸を時間刻み幅𝜏、縦軸を誤差推定∆𝑅であるグラフにすると図 36 になる。図 3-6 の黒い点線は𝑂(𝜏 )に従う直線である。

図3-6

1 次元 TDGL 方程式の誤差の時間刻み幅依存性(AFI-2)

図 3-6 を見ると AFI-2 を使用した 1 次元 TDGL 方程式の誤差の時間刻み幅依存性は𝑂(𝜏 )に従うことが分かる。これは図 3-1、図 3-4 から AFI-2 の誤差の時間刻み幅依存性は𝑂(𝜏 )に従っているように見える。ゆえに、AFI-2-opt の誤差の時間刻み幅依存性は𝑂(𝜏 )、AFI4-opt の誤差の時間刻み幅依存性は𝑂(𝜏 )と予想する。
AFI-2-opt と AFI-4-opt を使用した 1 次元 TDGL 方程式の誤差の時間依存性と誤差の時間刻み幅依存性を図に示す。AFI-2-opt を使用した 1 次元 TDGL 方程式の誤差の時間依存性を図 3-7、AFI-2-opt を使用した 1 次元 TDGL 方程式の誤差の時間刻み幅依存性を図 38、AFI-4-opt を使用した 1 次元 TDGL 方程式の誤差の時間依存性を図 3-9、AFI-4-opt を使用した 1 次元 TDGL 方程式の誤差の時間刻み幅依存性を図 3-10 に示す。図 3-8、図 310 の黒い点線は𝑂(𝜏 )に従う直線、図 3-10 の水色の点線は𝑂(𝜏 )に従う直線である。

38

図3-7

1 次元 TDGL 方程式の誤差の時間依存性(AFI-2-opt)

図3-8

1 次元 TDGL 方程式の誤差の時間刻み幅依存性(AFI-2-opt)

39

図3-9

1 次元 TDGL 方程式の誤差の時間依存性(AFI-4-opt)

図3-10

1 次元 TDGL 方程式の誤差の時間刻み幅依存性(AFI-4-opt)40

予想と同様に AFI-2-opt を使用した 1 次元 TDGL 方程式の誤差の時間刻み幅依存性は𝑂(𝜏 )に従った。しかし、予想とは異なり AFI-4-opt を使用した 1 次元 TDGL 方程式の誤差の時間刻み幅依存性は𝑂(𝜏 )ではなく𝑂(𝜏 )に従った。しかし、図 3-10 を見ると𝑡 = 0.2のグラフの𝜏 = 0.02から𝜏 = 0.01の間では𝑂(𝜏 )に従っている。
これは実際は AFI-2-opt の誤差の時間刻み幅依存性は𝑂(𝜏 )、AFI-4-opt の誤差の時間刻み幅依存性は𝑂(𝜏 )に従うが AFI-4-opt では初期条件から安定するまで変化が早くて捉えられず、時間刻み幅が変化したときに生まれる誤差に従っているからである。
これをはっきりさせるには、磁場や電流を印加し、オーダーパラメータ𝜓を変化させる必要があるだろう。

3.6 1 次元 TDGL 方程式の安定性
𝑚 = 1, 𝐿 = 10, 𝑁 = 100, ℎ = 0.1, 𝛾 = 1, 𝐴 = 0における時間刻み幅𝜏と𝛽(= −𝛼)を変更した
ときの周期的境界条件の 1 次元 TDGL 方程式の AFI の安定性を図 3-11 に示す。このとき
AFI-2、AFI-2-opt、AFI-4-opt を使用して数値計算を行った。横軸は時間刻み幅𝜏、縦軸は非線形パラメータ𝛽である。非線形パラメータ𝛽を設定したとき、丸が数値的安定性のある最大の時間刻み幅𝜏、クロスが数値的安定性のない最小の時間刻み幅𝜏である。

図3-11

1 次元 TDGL 方程式の AFI の安定性

まず、AFI-4-opt では𝛽が 0 から 170 の間では数値的安定性の境界線は時間刻み幅が 0.05
41

から 0.1 の間にあることが分かる。また𝛽が 180 から 200 の間では時間刻み幅が 0.05 以下
でか数値的に安定せず、𝛽が 200 を超えると時間刻み幅が 0.001 で数値シミュレーションを行っても不安定領域であった。
次に、AFI-2 について考える。AFI-2 では𝛽が 0 から 160 までは時間刻み幅が 0.2 以下では安定領域しか存在せず、𝛽が 170 から 200 の間で安定領域と不安定領域の境界線があることが分かった。
AFI-2-opt のときも AFI-2 と似たように𝛽が 0 から 110 までは時間刻み幅が 0.2 以下では安定領域しか存在せず、𝛽が 120 から 200 の間で安定領域と不安定領域の境界線があることが分かる。
これから、AFI それぞれを用いた時の 1 次元 TDGL 方程式の安定性にかなり違いがみられることが分かった。

42

第 5 章 結論
本研究では、ゲージ場下の量子現象を記述する時間発展方程式を高精度で解くことができる AFI の誤差評価を行った。
周期的境界条件を用いた 1 次元拡散方程式の誤差評価を行った。2 次のオーダーの LTS
分解を使用した AFI や 2 次のオーダーの最適化 LTS 分解を使用した AFI は、Euler 法やRunge-Kutta 法といった方法に比べて大きな時間刻み幅でも安定して数値計算が行えるが誤差は大きい。6 次のオーダーの最適化 LTS 分解を使用した AFI は、Euler 法や RungeKutta 法と比べても誤差の大きさは変わらないが、安定して数値計算行える時間刻み幅の大きさが変わらない。しかし、4 次のオーダーの最適化 LTS 分解を使用した AFI は、Euler 法や Runge-Kutta 法と比べても大きな時間刻み幅でも安定して数値計算が行え、誤差も Euler
法や Runge-Kutta 法と変わらず全体的に優れていた。
周期的境界条件を用いた 2 次元拡散方程式の誤差評価を行った。2 次のオーダーの LTS
分解を使用した AFI や 2 次のオーダーの最適化 LTS 分解を使用した AFI は、Euler 法やRunge-Kutta 法といった方法に比べて大きな時間刻み幅でも安定して数値計算が行えるが誤差は大きい。4 次のオーダーの最適化 LTS 分解を使用した AFI は、Euler 法や RungeKutta 法と比べても大きな時間刻み幅で安定して数値計算が行え、誤差も Euler 法や RungeKutta 法と変わらず全体的に優れていた。
ディリクレ境界条件を用いた 1 次元量子拡散方程式の誤差評価を行った。ベクトルポテ
ンシャルがないとき、2 次のオーダーの LTS 分解を使用した AFI や 2 次のオーダーの最適化 LTS 分解を使用した AFI は、Euler 法や Runge-Kutta 法といった方法に比べて大きな時
間刻み幅でも安定して数値計算が行えるが誤差は大きい。4 次のオーダーの最適化 LTS 分解を使用した AFI は、Euler 法や Runge-Kutta 法と比べても大きな時間刻み幅で安定して数値計算が行え、誤差も Euler 法や Runge-Kutta 法と変わらなかった。ベクトルポテンシ
ャルがあるとき、数値シミュレーションを実行した時、境界付近で解析解と数値解の振る舞
いが違くなり、Euler 法、Runge-Kutta 法や AFI といった方法に関係なく誤差の大きさは変
わらず、誤差の大きさは大きくなった。しかし、ベクトルポテンシャルに関わらず数値解法
それぞれの安定して数値計算できる時間刻み幅の大きさは変わらない。
周期的境界条件を用いた 1 次元量子拡散方程式の誤差評価を行った。時間刻み幅が小さいとき、印加するベクトルポテンシャルを大きくしていくのに従って誤差は大きくなっているが、ディリクレ境界条件を用いた 1 次元量子拡散方程式の誤差評価のときの誤差よりも十分に小さい。時間刻み幅が大きいときは誤差の大きさに違いは見られなかった。2 次のオーダーの LTS 分解を使用した AFI や 2 次のオーダーの最適化 LTS 分解を使用した AFIは、Euler 法や Runge-Kutta 法といった方法に比べて大きな時間刻み幅でも安定して数値計算が行えるが誤差は大きい。4 次のオーダーの最適化 LTS 分解を使用した AFI は、Euler 法
43や Runge-Kutta 法と比べても大きな時間刻み幅で安定して数値計算が行え、誤差も Euler 法や Runge-Kutta 法と変わらず全体的に優れていた。また、数値シミュレーションを実行した時の数値解と解析解の振る舞いから境界条件はディリクレ境界条件よりも周期的境界条件の方が誤差評価を行う際に適切だと言える。
周期的境界条件を用いた 1 次元時間依存 Ginzburg-Landau(TDGL)方程式の誤差評価を行った。2 次のオーダーの LTS 分解を使用した AFI を用いた 1 次元 TDGL の誤差の時間刻み幅依存性は 2 次のオーダーに従った。
4 次のオーダーの最適化 LTS 分解を使用した AFIを用いた 1 次元 TDGL の誤差の時間刻み幅依存性は 4 次のオーダーに従う。しかし、計算しているオーダーパラメータの値が一定の値で安定してしまうため、誤差推定に用いる 2 つの時間刻み幅で数値シミュレーションを実行した時に生まれる誤差が影響を与え、2 次のオーダーに従うように見えてしまう。つまり、1 次元 TDGL 方程式に AFI を用いて数値計算行った時の誤差は AFI に用いる LTS 分解のオーダーに従うことが分かる。
周期的境界条件を用いた 1 次元 TDGL 方程式の安定性を計算した。時間刻み幅と TDGL
方程式の展開係数𝛽の値を変更した時、数値的安定性があるかどうかで計算を行っている。
2 次のオーダーの LTS 分解を使用した AFI は𝛽が 0 から 160 の間では時間刻み幅が 0.2 以
下では不安定領域が存在せず、𝛽が 170 から 200 の間で不安定領域が存在しはじめ、安定領域と不安定領域の境界線があることがわかる。2 次のオーダーの最適化 LTS 分解を使用した AFI は𝛽が 0 から 110 の間では時間刻み幅が 0.2 以下では不安定領域が存在せず、𝛽が120 から 200 の間で不安定領域が存在しはじめ、安定領域と不安定領域の境界線があることがわかる。4 次のオーダーの最適化 LTS 分解を使用した AFI は𝛽が 0 から 170 の間では
時間刻み幅が 0.05 から 0.1 の間に安定領域と不安定領域の境界線があることが分かった。
また、𝛽が 180 から 200 の間では時間刻み幅が 0 から 0.05 の間に安定領域と不安定領域の境界線があることが分かる。これから、AFI それぞれを用いた時の 1 次元 TDGL 方程式の安定性にかなり違いがみられることが分かった。
AFI は Euler 法や Runge-Kutta 法といった既存の方法と比較して、拡散方程式、量子拡散方程式、ゲージ場存在下の TDGL 方程式の数値計算に用いても誤差の大きさはあまり違いは見られない。だが、Euler 法や Runge-Kutta 法では発散してしまうような大きな時間刻み
幅でも安定して数値計算を行うことができる。これは⾧時間必要な数値計算の回数が少な
くなり、数値シミュレーションの時間が短く済み、超伝導体内を可視化するうえでかなりの利点がある。
今後の課題として、4 次のオーダーの最適化 LTS 分解を使用した AFI に比べて 6 次のオーダーの最適化 LTS 分解を使用した AFI は小さな時間刻み幅でしか安定して数値計算が行えないため、大きな時間刻み幅でも安定して数値計算が行える高次のオーダーの LTS 分解を見つける必要がある。

44

謝辞
本研究に関して、様々な方々から多大なご助力をいただきました。
指導教員である小田部荘司先生には、研究に関する知識を何度もご教授頂き、また研究をするうえで、アドバイスや参考になる文献を紹介していただき、大変感謝しています。他にも海外の大学生や教授との学術的な交流の機会を下さり、研究以外にも貴重な体験をさせていただき大変感謝しています。
有明工業高等専門学校の松野哲也先生からは、何度も AFI や数値解析に関して指導をしていただきました。松野先生からいただいた AFI や数値解析に関する論文や資料は研究への理解や発展に役立ちました。誠に感謝いたします。
また、小田部研究室のメンバーには公私ともに助言をいただき、大変世話になりました。
この場をお借りして深くお礼申し上げます。

45

研究業績
1. 川畑 唯一、閻 洪、小田部 荘司、松野 哲也、馬渡 康徳 「2020 年応用物理学会九州支部学術 講演会」 2020 年 11 月 29 日

46

参考文献
[1] T. Matsuno, E.S.Otabe, and Y. Mawatari,” Explicit Integrators Based on a Bipartite
Lattice and a Pair of Affine Transformations to Solve Quantum Equations with Gauge Fields”,
J. Phys. Soc. Jpn. 89, 054006(2020).
[2] 松野哲也 著 量子化磁束動力学シミュレーション研究グループ

夏のセミナー 2020

資料 Ver.3
http://annex.jsap.or.jp/fluxoid/img/eqs_summer_v3.pdf(参照 2021 年 2 月 15 日)[3] T. Barthel and Y. Zhang, “Optimized Lie-Trotter-Suzuki decompositions for two andthree non-commuting terms”, Annals of Physics 418 (2020) 168165.
[4] 中田真秀 論文紹介：Lie-Trotter-鈴木分解について
https://www.slideshare.net/NakataMaho/lietrottersuzuki-159710114(参照 2021 年 2 月 15 日)[5] 伊理正夫・藤野和建 著「数値計算の常識」 共立出版社, 1985
[6] 松野哲也 小田部荘司 馬渡康徳 超電導体における電磁減少のシミュレーションのための陽的数値積分法 日本物理学会 2020 年秋季大会[7] 松野哲也 有明工業高等専門学校想像工学科：私信による(2020 年)[8] V. L. Ginzburg, L. D. Landau, Zh. Eksp. Teor. Fiz. 20 (1950) 1064
[9] 田中宏樹 アフィン変換数値積分法を用いた超伝導体内の量子化磁束線の可視化 九州
工業大学情報工学部学士学位論文(2020)[10] 濱田雄成 時間依存 Ginzburg-Landau 方程式を用いた様々な条件下における超伝導体の量子化磁束線の運動解析 九州工業大学情報工学部学士学位論文(2018) A-19

47