力定数の総和則をとりいれた最小二乗法の定式化 33

33

第 4 章 力定数の総和則をとりいれた最小

∑

i,t

2{(f(k_i, t;p₁,〜, pn)−y_i,t}∂f(k_i, t;p₁,〜, pn)

∂p_j′ = 0 (4.1.3)

である。ここでf(ki, t;p1,〜, pn)をあるforce constantの値で一次の項まで展開すると,

f(k_i, t;p₁,〜, pn) = f(k_i, t;p⁽⁰⁾₁ ,〜, p⁽⁰⁾_n ) +∑

j

∂f(ki, t;p⁽⁰⁾₁ ,〜, p⁽⁰⁾n )

∂p_j

¯¯¯¯

¯pj=p⁽⁰⁾_j

(p_j−p⁽⁰⁾_j ) (4.1.4) となる。ここで式4.1.3に式4.1.4を代入すると

∑

i,t,j

(∂f(k_i, t;p₁,〜, pn)

∂p_j

) (∂f(k_i, t;p₁,〜, pn)

∂p_j′

)

(p_j −p⁽⁰⁾_j ) (4.1.5)

=∑

i,t

(y_i,t−f(k_i, t;p₁,〜, pn))

(∂f(k_i, t;p₁,〜, pn)

∂p_j′

)

(j^′ = 1〜n) (4.1.6) となる。ここで四近接までの原子からの力の寄与を考えると一近接あたりϕ_r, ϕ_ti, ϕ_toと3 つのパラメータがあるから全部で3×4＝12のパラメータがある。また∆p_j =p_j −p⁽⁰⁾_j とし上式を行列の形にすると,













∑ i

(∂f(k_i, t;p₁,〜, p₁₂)

∂p₁

) (∂f(k_i, t;p₁,〜, p₁₂)

∂p₁ )

· · · · · · ∑

i

(∂f(k_i, t;p₁,〜, p₁₂)

∂p₁

) (∂f(k_i, t;p₁,〜, pn)

∂p₁₂ )

∑ i

(∂f(k_i, t;p₁,〜, p₁₂)

∂p₂

) (∂f(k_i, t;p₁,〜, p₁₂)

∂p₁ )

· · · · · · ∑

i

(∂f(k_i, t;p₁,〜, p₁₂)

∂p₂

) (∂f(k_i, t;p₁,〜, p_n)

∂p₁₂ )

.. .

.. .

.. .

.. . ..

.

.. .

.. .

.. . ..

.

.. .

.. .

.. .

∑ i

(∂f(k_i, t;p₁,〜, p_n)

∂p₁₂

) (∂f(k_i, t;p₁,〜, p_n)

∂p₁ )

· · · · · · ∑

i

(∂f(k_i, t;p₁,〜, p_n)

∂p₁₂ )(∂f(k_i, t;p₁,〜, p_n)

∂p₁₂ )





















∆p₁

∆p₂ .. . .. . .. .

∆p₁₂









=













∑ i

{y_i,t−f(x_i:p₁〜p₁₂)}(

∂f(k_i, t;p₁,〜, p_n)

∂p₁ )

∑ i

{y_i,t−f(x_i:p1〜p12)}(

∂f(k_i, t;p₁,〜, p_n)

∂p₂ )

.. . .. . .. .

∑ i

{y_i,t−f(x_i:p₁〜p₁₂)}(

∂f(k_i, t;p₁,〜, p_n)

∂p₁₂ )













(4.1.7)

ここで連立方程式4.1.7を解くことによって∆pj =p⁽¹⁾_j −p⁽⁰⁾_j の値を求めることができる。

p⁽ⁿ⁾_j はj番目のforce constantパラメータでn回目の対角化で求められたforce constantパ ラメータとしている。連立方程式を解くことで求めることができたp⁽¹⁾_j を初期値として式

第4章 力定数の総和則をとりいれた最小二乗法の定式化 35 (4.1.7)を繰り返し解くことによってp⁽²⁾_j を求めることができる。式4.1.1で定義をした実 験と計算の誤差の二乗の和Mが十分に小さくなるまで対角化を繰り返す。次に,3章で導 入した力定数の総和則のとりいれかたについて,第4近接までの原子の力の寄与を考えた 場合を例に説明する。はじめ,二種類の力定数の総和則,式3.2.9と式3.2.12を拘束条件と してラグランジュの未定乗数法をつかって試みたが式4.1.1の値を収束させることができ なかったため,式3.2.9と式3.2.12を

ϕ⁽⁴⁾_ti = − 1 14

(

ϕ⁽¹⁾_ti + 6ϕ⁽²⁾_ti + 4ϕ⁽³⁾_ti )

(4.1.8) ϕ⁽⁴⁾_to = − 1

14 (

ϕ⁽¹⁾_to + 6ϕ⁽²⁾_to + 4ϕ⁽³⁾_to )

(4.1.9) とし,一番遠くの原子で最も力の寄与が小さいと思われるパラメータの値を近い原子に関 するパラメータの値から決定されるようにした。

4.2 数値計算の手法

この節では実際にもちいて計算した手法やfittingをする際に必要なパラメータの選び 方について示す。

4.2.1 初期条件

前節でも説明したが最初に初期条件p⁽⁰⁾_j (j = 1,・・・n)を決め次に新しいパラメータ p¹_j (j = 1,・・・n)を得る。これを繰り返すことによって新しいパラメータを得ることがで きる。本研究ではR.Saitoらがグラフェンのフォノンの計算でもちいたforce constanrtの

値図4.1[13]をパラメータの値を最初の初期条件として用いた。

近接数 ϕ_r ϕ_ti ϕ_to  1  36.5 24.5 9.82  2  8.8 -3.23 -0.4  3  3.0 -5.25 0.15  4  -1.92 2.29 -0.58

図4.1: R.Saitoらがグラフェンのフォノンの計算でもちいたforce constanrtの値.

パラメータの初期値は実際に計算する際に非常に重要な要素であり,初期値によっては 収束しない場合もある。また近接数を5近接まで増やすときは第4近接で収束させた結果 を初期値として計算を繰り返していき同様の方法で14近接まで計算した。近接数を増や

4.2.2 _微係数

式4.1.7中の行列要素の関数f(k_i, t;p₁,〜, pn)をforce constantパラメータで微分する部 分∑

i

 

(∂f(k_i, t;p₁,〜, pn)

∂p_j

) (∂f(k_i, t;p₁,〜, pn)

∂p_j′

)

が含まれている。数値微分は前進差 分で次のように計算できる。

∂f(k_i, t;p_j)

∂p_j = f(k_i, t;p_j +δp_j)−f(k_i, t;p_j)

δp_j (4.2.1)

しかし経験的に前進差分で微分を扱うと誤差が大きくなってしまうため実際の計算ではに 5点法で計算した。

∂f(k_i, t;p_j)

∂p_j = f(x_i, p_j−2δp_j)−8f(x_i, p_j−δp_j) + 8f(x_i, p_j +δp_j)−f(x_i, p_j + 2δp_j) 12δp_j

(4.2.2) ここでδp_jの値はパラメータの種類によって非常に敏感になるものがあり選び方によって は発散してしまうものもあるため,δpj の値を変化させ収束のしかたを調べる必要がある。

4.2.3 新しいパラメータの決定

定式化の節でも書いたように行列を対角化することによってδp_j が求まり

p^new_j =p^old_j +δp_j (4.2.3)

によって新しいパラメータを求めることができるが,定式化のところで定義をしたMが収 束するまで非常に時間がかかることがあったため次に記すように改良した。実数λをもち いて新しいパラメータp^new_j を以下のようにした。

p^new_j =p^old_j +λδp_j (4.2.4) ここでλの値を0.1刻みで10まで変化させ実験と計算との誤差の値M =∑

i

{f(ki, t;p1〜pn)−yi,t}² が最小になるようなλの値を用いることにした。このことによりMを収束させるまでの

計算時間を短縮させることができる。

37

第 5 _{章 計算結果} I( _{力定数の総和則をとり}

ドキュメント内 ナノグラファイトの振動構造 (ページ 34-38)