21.1 1 次元の場合

例題1

√2，√ 3，√

4，…，√

10を順に求めるプログラムを作成する．具体的には，y=x² として，y= 2,3,…,10を満たすx >0を順次求めていく．

実は，幾何学的非線形性を厳密に考慮した仮想仕事式を有限要素離散化して得ら れる．非線形の連立方程式（以下では平衡方程式と呼ぶ）も，外力などの境界条 件によって定まるaに対して，平衡方程式f(x)が

230 第21章 非線形解法の代表 Newton-Raphson法

f(x) =a (21.1)

を満たすようなxを求めるという形式の問題である．このaに対する解を直接求 めることができる場合もあるが，f(x)の非線形性が強いと，aに対する解を求め るために，たとえば，aを10等分して，0< a₁ < a₂ <… < a₁₀ =aとして，まず はa₁に対する解x₁を求め，次にこれをもとに，a₂に対する解x₂を求め，順にx₉ をもとにa₁₀(=a)に対する解xを求める．このように与えられたaに対して中間 的なa_iに対して順に解を求めるのが増分解析の基本である．

1.1収束判定

実際に解を求めるプロセスの説明をする前に，近似解の意味について考えてみ る．まず，近似的であれ何であれ，解の候補x⁽ⁱ⁾が，解としてふさわしいかどうか 判定する必要がある．このために，下式に基づき，残差rを定義する．

r =a−f! x⁽ⁱ⁾"

(21.2) 当然ながら，正解xは，r= 0となる．近似解x⁽ⁱ⁾では，x⁽ⁱ⁾ xであれば，|r| 1 となるはずである．ここで注意すべきは，計算機で計算する限り，r= 0となる保 証がないということである．例えば，a=√

2と，値を代入したところで，aに無 理数としての√

2を代入できるわけではない．あくまで，√

2を有限桁数の浮動小 数点で近似したものが代入されているだけで，OSやコンパイラが変われば，その 値自体も異なる場合もある．そのため，何らかのアルゴリズムで，いかに正確に

√2の近似解を求めたとしても，計算機の内部表現の意味での√

2に，bit単位で一 致し，r= 0となることは保証されないのである．そのため通常は，を十分小さ い正数(10⁻⁸〜10⁻¹⁰など)として，|r|< が満たされれたら，x⁽ⁱ⁾はxを十分な精 度で近似しているとみなす．当然ながら，は小さいほど近似の精度が上がること が期待できるが，倍精度の実数でも約15桁の有効数字しかないので， = 10⁻¹⁰⁰ としも意味が無い．また， = 10⁻⁶〜10⁻⁸ 程度でも，同じ基準で判定しても後続 の計算で著しく反復回数が増加するなどの困難を生じることなく，解き進むこと ができるのであれば，実質的に十分な精度の解を求められていると考えられる場 合は多い．このように，ある与えられたに対して，|r| < を満たすかどうか判 定することを収束判定と呼び，得られた近似解x⁽ⁱ⁾を収束解と呼ぶ．

1.2モンテカルロ法

残差が十分小さい近似解x⁽ⁱ⁾を収束解とするのであれば、最も単純なアルゴリ ズムとしてΔxを経験などを元に設定して

x⁽ⁱ⁾_k+1 =x_k+ Δx・i (21.3)

21.1. 1次元の場合 231 のようにk段階目の収束解x_kを元にk+ 1段階目の収束解x_k+1を求めるというア ルゴリズムを考えることができる。あるいはa⁽ⁱ⁾を適当な乱数として、

x⁽ⁱ⁾_k+1 =x_k+ Δx・a⁽ⁱ⁾ (21.4) としてもよい。このような方法であっても、反復回数を十分大きくとればそれな りの近似解を得ることができる。後者の乱数を用いるアルゴリズムは特にモンテ カルロ法と呼ばれている。後述のNewton-Raphson法と比較して長所は、どのよ うなf(x)であっても、即ちたとえ連続な関数でなく、ロジックによって離散的に f(x)の値が得られるような場合でもf(x)さえ計算できれば収束解を得ることがで きることにある。短所は、もちろん実際に実行してみるとわかるように一般に非 常に時間がかかることである。

2.Newton-Raphson法

もしf(x)を求めることができればモンテカルロ法よりはるかに効率よく解を 求める方法として、Newton-Raphson法が知られており、非線形有限要素法でも

Newton-Raphson法や、これを元にした修正Newton法や、割線法がよく用いられ

ている。 

微分可能な関数は微小区間[x, x+dx]の間では線形とみなすことができるという のが微分の基本であるが、有限な、即ち微小ではないΔxについてもある程度小 さければ区間[x, x+ Δx]においてf(x)はそれほど線形から外れることはない。そ こで与えられたa₁に対するf(x) = a₁ の近似解x₁ が既に求められているとき、

f(x) =a₂に対する近似解x₂の題１次予測としてx⁽¹⁾₂ を以下のように求める。

f(x₁)Δx⁽¹⁾₁ = a₂−f(x₁) (21.5)

→Δx⁽¹⁾₁ = 1

f(x₁)(a₂−f(x₁)) (21.6) x⁽¹⁾₂ =x₁+ Δx⁽¹⁾₁ (21.7)

232 第21章 非線形解法の代表 Newton-Raphson法

a

₁

a

₁

a

₂

a

₂

x

₁

x

₁

x

⁽¹⁾₂

x

⁽²⁾₂

x

⁽¹⁾₂

Δ x

⁽¹⁾₁

Δ x

⁽²⁾₂

f ( x

⁽¹⁾₂

) f ( x

⁽¹⁾₂

) f ( x

⁽²⁾₂

)

a

₂

− f ( x

₁

)

a

₂

− f ( x

⁽¹⁾₂

)

このときの残差r⁽¹⁾₂ は以下のように求められる。

r₂⁽¹⁾=a₂−f(x⁽¹⁾₂ ) (21.8) もし、r₂⁽¹⁾が十分に小さければこれで計算を終えることも不可能ではないが、多く の場合は、十分に解が収束している。即ち、これ以降計算を行ったときにrが減 少するかどうか定かでないので、次式に従い第２段階での近似x⁽²⁾₂ を求める。

f(x⁽¹⁾₂ )Δx⁽²⁾₂ = a₂−f(x⁽¹⁾₂ ) (21.9)

→Δx⁽²⁾₂ = 1 f(x⁽¹⁾₂ )

(a₂−f(x⁽¹⁾₂ )) (21.10)

x⁽²⁾₂ =x⁽¹⁾₂ + Δx⁽²⁾₂ (21.11) このプロセスを図示する。これ以降順に以下のように近似解x⁽ⁱ⁾₂ を元にx⁽ⁱ⁺¹⁾₂ を求 める。

f(x⁽ⁱ⁾₂ )Δx⁽ⁱ⁺¹⁾₂ = a₂−f(x⁽ⁱ⁾₂ ) (21.12)

→Δx⁽ⁱ⁺¹⁾₂ = 1

f(x⁽ⁱ⁾₂ )(a₂−f(x⁽ⁱ⁾₂ )) (21.13)

x⁽ⁱ⁺¹⁾₂ = x⁽ⁱ⁾₂ + Δx⁽ⁱ⁺¹⁾₂ (21.14)

r⁽ⁱ⁺¹⁾₂ = a₂−f(x⁽ⁱ⁺¹⁾₂ ) (21.15)

|r₂⁽ⁱ⁺¹⁾|< εとなれば反復計算を終了することになる。

21.1. 1次元の場合 233 以下に説明したNewton-Raphson法により√

2の値を手計算で求めてみよう。ま ずy=x²を満たす初期値として、x₁ = 1, y₁ = 1からスタートすることにしよう。

実は、x= 0, y = 0からスタートすることが多い。これは平衡方程式においては、

外力が0なら無変形ということに対応する自明な解だからである。たまたま、今 回の問題はf(0) = 0となり、Newton-Raphson法の最初で求積できないので、例 外的にx = 0の初期値を採用している。ただし、実際の物理現象でも何らかの都

合でf(0) = 0となることがあり、この種の工夫が必要になることは少なくない。

f(x) =x², df(x)

dx = 2x (21.16)

a₂ = 2, x₁ = 2を代入して、

df(1)

dx = 2, 2Δx⁽¹⁾₁ = 2−1 (21.17)

→Δx⁽¹⁾₁ = 1

2 (21.18)

x⁽¹⁾₂ =x₁+ Δx⁽¹⁾₁ = 1 +1 2 = 3

2 = 1.5 (21.19)

a₂ = 2, x⁽¹⁾₂ = ³₂ を代入して、

df(³₂)

dx = 3, f(3 2) = 9

4, 3Δx⁽²⁾₂ = 2− 9

4 (21.20)

→Δx⁽²⁾₂ =− 1

12 (21.21)

x⁽²⁾₂ =x⁽¹⁾₂ + Δx⁽²⁾₂ = 3 2 − 1

12 = 17

12 = 1.4166· · · (21.22) a₂ = 2, x⁽²⁾₂ = ¹⁷₁₂を代入して、

df(¹⁷₁₂) dx = 17

6 , f(17

12) = 289

144, 17

6 Δx⁽³⁾₂ = 2− 289

144 (21.23)

→Δx⁽³⁾₂ =− 1

408 (21.24)

x⁽³⁾₂ =x⁽²⁾₂ + Δx⁽³⁾₂ = 17 12− 1

408 = 577

408 = 1.4142156· · · (21.25) a₂ = 2, x⁽³⁾₂ = ⁵⁷⁷₄₀₈を代入して、

df(⁵⁷⁷₄₀₈)

dx = 577

204, f(577

408) = 577·577

408·408, 577

204Δx⁽⁴⁾₂ = 2− 577·577

408·408 (21.26)

234 第21章 非線形解法の代表 Newton-Raphson法

→Δx⁽⁴⁾₂ =− 1

470832 (21.27)

x⁽⁴⁾₂ =x⁽³⁾₂ + Δx⁽⁴⁾₂ = 577

408 − 1

470832 = 665857

470832 = 1.41421356237· · · (21.28) 反復数4回でも精度の高い近似解が得られていることが分かる。

このプロセスをループで表してみよう。

dok = 2, 10 a_k =k

doi = 1, 100

df_k⁽ⁱ⁻¹⁾

dx Δx⁽ⁱ⁾_k =a_k−f(x⁽ⁱ⁻¹⁾_k )→Δx⁽ⁱ⁾_k x⁽ⁱ⁾_k =x⁽ⁱ⁻¹⁾_k + Δx⁽ⁱ⁾_k

r⁽ⁱ⁾_k =a_k−f(x⁽ⁱ⁾_k ) if |r⁽ⁱ⁾_k |< ε

convergence flag = 1 exit loop

else

convergence flag = 0 end if

end

if convergence flag = 1 write x

else

write ” not converged ” exit

end if end

実際のプログラムではループカウンターにi、kなどよく使われる文字を使うと途 中で紛れてbugになる可能性があるので、kの代わりにiincre(index of increment)、

iの代わりにiiter(index of interation)を用いることにする。また、ループの最大 値もこれにあわせ、nincre、niterとすれば、データファイルからの読み込みに対 応できる。

このサンプルプログラムを示す。

1dim.f ところで、f(x)が一次関数、たとえばf(x) =cx(cは定数)であったら、

ドキュメント内 all.dvi (ページ 44-50)