これまでのまとめ ( 学年末試験に向けて )

これまでのまとめ ( 学年末試験に向けて )

山本昌志

^∗ 2005

2

14

1 はじめに

学年末試験に向けて、学習のポイントを示す。後期の中間試験から、これまで以下について数値計算の学 習を行った。

•

(反復法) –

–

– SOR

•

–

–

•

–

–

–

•

–

–

ここでは、これらについて、簡単にまとめてある。これが学習の最重要ポイントであるのでしっかり勉強し てほしい。このプリントを見ながら、分からない部分は授業中配ったプリントで補い、理解することが学習 のこつである。

数値計算のいろいろな方法を学習したが 、ほとんど 全ての方法は

1

それを使いこなせれば 、ここでの授業での学習は完璧である。学年末試験では、使いこなせるか否かを調べ ることは難しいので、式を導くことを中心に出題する。

∗

2 連立一次方程式 ( 反復法 )

実用的なプログラムでは、非常に大きな連立方程式を計算しなくてはならない。数百万元に及ぶことも珍 しくない。これを、ガウス・ジョルダン法で計算するの時間的にほとんど 不可能である。そこで、これより は格段に計算の速い方法が用いられる。ここでは、その一つとして反復法を簡単に説明する。

当然ここでも、連立方程式

Ax = b (1)

を満たす

x

x

ここで、ある計算により

n

x n

n→∞ lim x n = x (2)

になったとする。この様に計算回数を増やして、真の解に近づける方法を反復法という。

この様な方法は、行列

A

S − T

Sx k+1 = T x k + b (3)

とすればよい。ここで、

x k

α

(3)

(1)

x

(3)

(1)

_k

x

2.1

計数行列

A

S

(Jacobi)



 

 

 

 

a 11 a 12 a 13 . . . a 1n

a 21 a 22 a 23 . . . a 2n

a 31 a 32 a 33 . . . a 3n

.. . .. . .. . . .. .. . a n1 a n2 a n3 . . . a nn



 

 

 

 

=



 

 

 

 

a 11 0 0 . . . 0 0 a 22 0 . . . 0 0 0 a 33 . . . 0 .. . .. . .. . . .. .. . 0 0 0 . . . a nn



 

 

 

  +



 

 

 

 

0 a 12 a 13 . . . a 1n

a 21 0 a 23 . . . a 2n

a 31 a 32 0 . . . a 3n

.. . .. . .. . . .. .. . a n1 a n2 a n3 . . . 0



 

 

 

  (4)

と分解する。右辺第

1

S

2

−T

k+1

S

S ⁻¹ =



 

 

 

 

a ⁻¹ ₁₁ 0 0 . . . 0 0 a ⁻¹ ₂₂ 0 . . . 0 0 0 a ⁻¹ ₃₃ . . . 0 .. . .. . .. . . .. .. . 0 0 0 . . . a ⁻¹ _nn



 

 

 

 

(5)

と簡単である。k+1番目の近似解は、x

k+1 = S ⁻¹ (b + T x k )

k

x ^(k) ₁ , x ^(k) ₂ , x ^(k) ₃ , · · · , x ^(k) _n

x ^(k+1) ₁ = a ⁻¹ ₁₁ n

b 1 −

³

a 12 x ^(k) ₂ + a 13 x ^(k) ₃ + a 14 x ^(k) ₄ + · · · + a 1n x ^(k) _n

´o

x ^(k+1) ₂ = a ⁻¹ ₂₂ n b ₂ − ³

a ₂₁ x ^(k) ₁ + a ₂₃ x ^(k) ₃ + a ₂₄ x ^(k) ₄ + · · · + a _2n x ^(k) _n ´o x ^(k+1) ₃ = a ⁻¹ ₃₃ n

b 3 − ³

a 31 x ^(k) ₁ + a 32 x ^(k) ₂ + a 34 x ^(k) ₄ + · · · + a 3n x ^(k) _n ´o .. .

x ^(k+1) _n = a ⁻¹ _nn n b n − ³

a n1 x ^(k) ₁ + a n2 x ^(k) ₂ + a n3 x ^(k) ₃ + · · · + a nn−1 x ^(k) _n−1 ´o

(6)

と計算できる。これが 、ヤコビ法である。

2.2

ヤコビ法の特徴では、x

^(k+1)

x ^(k)

x ^(k+1)

x ^(k)

個人で大きなメモリーを使い計算することは許されるが 、ちょっと前まではできるだけメモリーを節約した プログラムを書かなくてはならなかった。

そこで、

x ^(k+1)

ち、x

^k+1 _i

x ^k+1 _i = a ⁻¹ _ii n

b i − (a i1 x ^(k+1) ₁ + a i2 x ^(k+1) ₂ + a i3 x ^(k+1) ₃ + · · · + a ii−1 x ^(k+1) _i−1 +

a ii+1 x ^(k) _i+1 + a ii+2 x ^(k) _i+2 + a ii+3 x ^(k) _i+3 + · · · + a in x ^(k) _n ) o

(7)

x ^(k+1) ₁ = a ⁻¹ ₁₁ n

b 1 −

³

a 12 x ^(k) ₂ + a 13 x ^(k) ₃ + a 14 x ^(k) ₄ + · · · + a 1n x ^(k) _n

´o

x ^(k+1) ₂ = a ⁻¹ ₂₂ n

b 2 −

³

a 21 x ^(k+1) ₁ + a 23 x ^(k) ₃ + a 24 x ^(k) ₄ + · · · + a 2n x ^(k) _n

´o

x ^(k+1) ₃ = a ⁻¹ ₃₃ n b ₃ − ³

a ₃₁ x ^(k+1) ₁ + a ₃₂ x ^(k+1) ₂ + a ₃₄ x ^(k) ₄ + · · · + a _3n x ^(k) _n ´o .. .

x ^(k+1) _n = a ⁻¹ _nn n

b n −

³

a n1 x ^(k+1) ₁ + a n2 x ^(k+1) ₂ + a n3 x ^(k+1) ₃ + · · · + a nn−1 x ^(k+1) _n−1

´o

(8)

と計算できる。これが 、ガウス・ザイデル法である。

2.3 SOR

ガウス・ザイデル法をもっと改善する方法がある。ガウス・ザイデル法の解の修正は、x

k+1 − x k

ω

k+1 − x k )

(SOR

Successive Over-Relaxation)

具体的な計算手順は、次のようにする。ここでは、ガウス・ザイデル法の式

(8)

x ˜ ^(k+1) _i

˜

x ^(k+1) ₁ = a ⁻¹ ₁₁ n b ₁ − ³

a ₁₂ x ^(k) ₂ + a ₁₃ x ^(k) ₃ + a ₁₄ x ^(k) ₄ + · · · + a _1n x ^(k) _n ´o x ^(k+1) ₁ = x ^(k) ₁ + ω ³

˜

x ^(k+1) ₁ − x ^(k) ₁ ´

˜

x ^(k+1) ₂ = a ⁻¹ ₂₂ n b 2 − ³

a 21 x ^(k+1) ₁ + a 23 x ^(k) ₃ + a 24 x ^(k) ₄ + · · · + a 2n x ^(k) _n ´o x ^(k+1) ₂ = x ^(k) ₂ + ω ³

˜

x ^(k+1) ₂ − x ^(k) ₂ ´

˜

x ^(k+1) ₃ = a ⁻¹ ₃₃ n

b 3 −

³

a 31 x ^(k+1) ₁ + a 32 x ^(k+1) ₂ + a 34 x ^(k) ₄ + · · · + a 3n x ^(k) _n

´o

x ^(k+1) ₃ = x ^(k) ₃ + ω

³

˜

x ^(k+1) ₃ − x ^(k) ₃

´

.. .

˜

x ^(k+1) _n = a ⁻¹ _nn n

b n −

³

a n1 x ^(k+1) ₁ + a n2 x ^(k+1) ₂ + a n3 x ^(k+1) ₃ + · · · + a nn−1 x ^(k+1) _n−1

´o

x ^(k+1) _n = x ^(k) _n + ω

³

˜

x ^(k+1) _n − x ^(k) _n

´

(9)

これが 、SOR法である。

ここで、問題なのが加速緩和係数

ω

1

ザイデル法で収束しないような問題には使える。

従って、1以上の値にしたいわけであるが 、余り大きくすると、発散するのは目に見えている。これにつ いては、2を越えると発散することが分かっている。最適値となると、だいたい

1.9

3 補間法

実験やシミュレーションを行うと、離散的にデータが得られるのは普通である。例えば 、半導体の電圧・

電圧特性の測定では、(V

0 , I 0 ), (V 1 , I 1 ), (V 2 , I 2 ), (V 3 , I 3 ), · · · , (V N , I N )

このデータはグラフ化して解析を進める。このデータの場合、2次元のグラフ上に測定点が並ぶことは、今 まで学習してきたとおりである。

実験等を通して得られる結果は離散的であるが、実際の現象は連続的なことが多い。この離散的な値を用 いて、測定点の間の値、ここでは電流と電圧の関係を求めるのが補間法の役割である。ここで学習したラグ ランジュ補間もスプライン補間も、全てのグラフ上の測定点を通る曲線の方程式を求めている。

2

(x, y)

(x, y)

3.1

平面座標上に

N + 1

N

¹

2

この性質を利用すると、N

+ 1

3

(x 0 , y 0 ), (x 1 , y 1 ), (x 2 , y 2 )

2

y = ax ² + by + c

a, b, c

x

y

a, b, c

コンピューターを用いて、N

+ 1

N

N + 1

N

補間という目的からすると、関数を形成する係数なんか、全く興味の対象外なのである。そこで、係数が分 からなくても、N次関数を示すものとして、ラグランジュ補間が使われる。

2

N + 1

0 , y 0 ) , (x 1 , y 1 ), (x 2 , y 2 ), · · · , (x N , y N )

y = (x − x 1 )(x − x 2 )(x − x 3 ) · · · (x − x N )

(x ₀ − x ₁ )(x ₀ − x ₂ )(x ₀ − x ₃ ) · · · (x ₀ − x _N ) y 0 + (x − x 0 )(x − x 2 )(x − x 3 ) · · · (x − x N ) (x ₁ − x ₀ )(x ₁ − x ₂ )(x ₁ − x ₃ ) · · · (x ₁ − x _N ) y 1

+ (x − x 0 )(x − x 1 )(x − x 3 ) · · · (x − x N )

(x 2 − x 0 )(x 2 − x 1 )(x 2 − x 3 ) · · · (x 2 − x N ) y 2 + (x − x 0 )(x − x 1 )(x − x 2 ) · · · (x − x N ) (x 3 − x 0 )(x 3 − x 1 )(x 3 − x 2 ) · · · (x 3 − x N ) y 3

· · · + (x − x 0 ) · · · (x − x k−1 )(x − x k+1 ) · · · (x − x N )

(x k − x 0 ) · · · (x k − x k−1 )(x k − x k+1 ) · · · (x k − x N ) y k + · · · + (x − x 0 )(x − x 1 )(x − x 2 ) · · · (x − x N −1 )

(x N − x 0 )(x N − x 1 )(x N − x 2 ) · · · (x N − x N−1 ) y N

(10)

この式

(10)

•

N

N

N

(N

• x

x 0 , x 1 , x 2 , · · · , y N

y

y 0 , y 1 , y 1 , · · · , x N

(x 0 , y 0 ) , (x 1 , y 1 ), (x 2 , y 2 ), · · · , (x N , y N )

となっている。

1 x

y

(10)

L(x) =

X N

k=0

L k (x)y k

ただし 、

L k (x) =

N (j6=k) Y

j=0

x − x j

x k − x j

(11)

となる。

ラグランジュ補間の考え方は単純で、その計算も簡単である。しかし 、補間の点数が増えてくると、ラグ ランジュの補間には問題が生じる。ラグランジュの補間では、補間の点数が増えてくると大きな振動が発生 して、もはや補間とは言えなくなる。ラグランジュの補間には常にこの問題が付きまうので、データ点数が 多い場合は使えなくなる。

3.2

3.2.1

ラグランジュの補間は、データ点数が増えてくると関数が振動し問題が発生する。そこで、補間する領域 をデータ間隔

[x i , x i+1 ]

(spline interpolation)

補間をするデータは 、先と同じ ように

(x 0 , y 0 ) , (x 1 , y 1 ), (x 2 , y 2 ), · · · , (x N , y N )

[x j , x j+1 ]

S j (x)

1

x _j x _j-1

x _j-2 x _j+1 x _j+2 x _j+3

S _j S _j-1 S _j-2

S _j+1 S _j+2

x y

図

1:

3.2.2

3

S j (x) = a j (x − x j ) ³ + b j (x − x j ) ² + c j (x − x j ) + d j (j = 0, 1, 2, 3, · · · , N − 1) (12)

a j , b j , c j , d j

これらの

4N

•

S j (x)

•

1

− 1

•

2

− 1

以上の条件を課すと

4N − 2

4N

この不足を補うために、いろいろな条件が考えられるが 、通常は両端

x 0

x N

2

0

₀ ⁰⁰ (x 0 ) = S _N ⁰⁰ ₋₁ (x N ) = 0

(natural spline)

1

これで全ての条件が決まった。あとは、この条件に満たす連立方程式を求めるだけである。

4 数値積分

4.1

定積分、

S = Z _b

a

f (x)dx (13)

の近似値を数値計算で求めることを考える。積分の計算は面積の計算であるから 、図

2

[a, b]

N

T

T = h f (a) + f (a + h)

2 + h f (a + h) + f (a + 2h)

2 + h f (a + 2h) + f(a + 3h)

2 + · · ·

+ h f (a + (N − 1)h) + f (a + N h) 2

= h 2

N X −1

j=0

[f (a + jh) + f(a + (j + 1)h)]

(14)

となる。これが数値積分の台形公式である。なんのことはない、積分を台形の面積に置き換えているだけで ある。

図

2:

4.2

台形公式の考え方は簡単であるが、精度はあまりよくない。そこで、よく似た考え方で精度が良いシンプ ソンの公式を説明する。台形公式は、分割点の値を一次関数

(直線)

2

(x j , x j+1 , x j+2 )

2

P(x)

P (x) = (x − x j+1 )(x − x j+2 )

(x j − x j+1 )(x j − x j+2 ) f (x j ) + (x − x j )(x − x j+2 )

(x j+1 − x j )(x j+1 − x j+2 ) f (x j+1 )

+ (x − x j )(x − x j+1 )

(x j+2 − x j )(x j+2 − x j+1 ) f (x j+2 ) (15)

3

これを、区間

[x j , x j+2 ]

(15)

図

3:

[x j , x j+2 ]

2

第

1

式

(15)

1

= Z _{x j+2}

x j

(x − x j+1 )(x − x j+2 )

(x j − x j+1 )(x j − x j+2 ) f (x j )dx

= Z _2h

0

(x j + ξ − x j+1 )(x j + ξ − x j+2 )

(x j − x j+1 )(x j − x j+2 ) f (x j )dξ

= Z _2h

0

(ξ − h)(ξ − 2h)

(−h)(−2h) f (x j )dξ

= f (x j ) 2h ²

Z _2h

0

ξ ² − 3hξ + 2h ² dξ

= f (x j ) 2h ²

· ξ ³

3 − 3h ξ ² 2 + 2h ² ξ

¸ 2h

0

= h 3 f(x _j )

(16)

同様に、第

2,3

式

(15)

2

= 4h

3 f (x j+1 ) (17)

式

(15)

3

= h

3 f (x j+2 ) (18)

となる。以上より、近似した

2

P(x)

[x j , x j+2 ]

Z _{x j+2}

x j

P(x)dx = h

3 {f (x j ) + 4f (x j+1 ) + f (x j+2 )} (19)

これは、ある区間

[x j , x j+2 ]

2h

[a, b]

S

(19)

を足し合わせればよい。ただし 、j

= 0, 2, 4, 6

S = h

3 {f (x 0 ) + 4f (x 1 ) + f (x 2 )} + h

3 {f (x 2 ) + 4f (x 3 ) + f (x 4 )} + h

3 {f (x 4 ) + 4f (x 5 ) + f (x 6 )} + · · ·

· · · + h

3 {f (x N −2 ) + 4f (x N−1 ) + f (x N )}

= h

3 {f (x 0 ) + 4f (x 1 ) + 2f (x 2 ) + 4f (x 3 ) + 2f (x 4 ) + 4f(x 5 ) + 2f (x 6 ) + · · ·

· · · + 2f (x N −2 ) + 4f (x N−1 ) + f (x N )}

(20)

N ⁴

この式から、分割数

N

4.3

積分の境界が複雑な場合、乱数を使うモンテカルロ積分が適している。例えば 、関数

f(x 1 , x 2 , x 3 , ·, x M )

Z

Ω

f (x 1 , x 2 , x 3 , · · · , x M )dx 1 dx 2 dx 3 · · · dx M = V Ω hf i (21)

V

M

Ω

i

3

(21)

積分を計算するために、まずは体積である。これは体積の計算が容易な単純な形状の内部に、領域

Ω

N

Ω

N Ω

V r

Ω

V Ω

V Ω ' N Ω

N V r (22)

の関係がある。右辺はコンピューターにより容易に計算できる。ランダムな点の数

N

残りは、体積内部の平均

hf i

Ω

hf i ' 1 N Ω

X

i

f (r _i )

Ω

(23)

r i

i

(x 1 , x 2 , x 3 , ·, x M )

Ω

以上より、モンテカルロ法を用いると、体積

V

hf i

(21)

5 差分法による偏微分方程式の数値計算

5.1

2

∂ ² φ

∂x ² + ∂ ² φ

∂y ² = 0 (24)

を数値計で解くことを考える。まずは、いつものように、解

φ(x, y)

y

±h

φ(x + h, y) = φ(x, y) + ∂φ

∂x h + 1 2!

∂ ² φ

∂x ² h ² + 1 3!

∂ ³ φ

∂x ³ h ³ + 1 4!

∂ ⁴ φ

∂x ⁴ h ⁴ + · · · (25) φ(x − h, y) = φ(x, y) − ∂φ

∂x h + 1 2!

∂ ² φ

∂x ² h ² − 1 3!

∂ ³ φ

∂x ³ h ³ + 1 4!

∂ ⁴ φ

∂x ⁴ h ⁴ − · · · (26)

∂ ² φ

∂x ²

¯ ¯

¯ ¯

x,y

= 1

h ² [φ(x + h, y) − 2φ(x, y) + φ(x − h, y)] − O(h ² ) (27)

h

²

(27)

O(h ² )

y

∂ ² φ

∂y ²

¯ ¯

¯ ¯

x,y

= 1

h ² [φ(x, y + h) − 2φ(x, y) + φ(x, y − h)] − O(h ² ) (28)

これらの式

(27)

(28)

2

(24)

φ(x + h, y) + φ(x − h, y) + φ(x, y + h) + φ(x, y − h) − 4φ(x, y) = 0 (29)

(x, y)

の値

φ(x, y)

実際にこの式を数値計算する場合、計算領域を間隔

h

²

y

h

(29)

y

(29)

実際、数値計算をする場合、φ(x, y)や

φ(x + h, y)

φ(x, y) = φ(ih, jh)

= U i j

(30)

(29)

U i+1 j + U i−1 j + U i j+1 + U i j−1 − 4U i j = 0 (31)

2

となり、数値計算し易い形になる。

ラプラス方程式は式

(31)

連立方程式を解くわけであるが、このままでは、式の数と未知数の数が異なる。格子点でのポテンシャル の値を求めるためには、境界条件を設定する必要がある。それにより、式の数と未知数の数が同一になり、

格子点でのポテンシャルを求めることができる。

5.2

弦の長さが

1、そこを伝わる波の速度を 1

 

 

 



 

 

 

∂ ² u

∂x ² = ∂ ² u

∂t ² (0 5 x 5 1, 0 5 t)

u(x, 0) = φ(x), ∂u

∂t (x, 0) = ψ(x) (0 5 x 5 1) u(0, t) = u(1, t) = 0

(32)

となる。弦を伝わる波の速度は

1、弦の長さも 1

波動方程式の他に、初期条件と境界条件がある。力学的状態は、ある時刻、ここでは

t = 0

4

弦

固 定 点

方 向 の速 度

図

4:

t = 0

(スナップショット)。初期条件と境界条件が表されており、y

ψ(x)

まずは、波動方程式を差分方程式に書き直すことからはじめる。これも、いつものように、解

u(x, t)

∆x、時間軸方向の微小変位を ∆t

u(x + ∆x, t) = u(x, t) + ∂u

∂x ∆x + 1 2!

∂ ² u

∂x ² (∆x) ² + 1 3!

∂ ³ u

∂x ³ (∆x) ³ + 1 4!

∂ ⁴ u

∂x ⁴ (∆x) ⁴ + · · · u(x − ∆x, t) = u(x, t) − ∂u

∂x ∆x + 1 2!

∂ ² u

∂x ² (∆x) ² − 1 3!

∂ ³ u

∂x ³ (∆x) ³ + 1 4!

∂ ⁴ u

∂x ⁴ (∆x) ⁴ − · · ·

(33)

となる。これらの式の辺々を足し合わせえると、

∂ ² u

∂x ²

¯ ¯

¯ ¯

x,t

= 1

∆x ² [u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t)] − O(∆x ² ) (34)

∆x

²

( 34)

1

∂ ² u

∂t ²

¯ ¯

¯ ¯

x,t

= 1

∆t ² [u(x, t + ∆t) − 2u(x, t) + u(x, t − ∆t)] − O(∆t ² ) (35)

これらの式

(34)

(35)

(32)

1

∆x ² [u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, t)] = 1

∆t ² [u(x, t + ∆t) − 2u(x, t) + u(x, t − ∆t)] (36)

u(x, t + ∆t) = 2u(x, t) − u(x, t − ∆t) + ∆t ²

∆x ² [u(x + ∆x, t) − 2u(x, t) + u(x − ∆x, y)] (37)

t

t − ∆t

t + ∆t

(37)

t

t − ∆t

+ ∆t

実際に式

(37)

∆x、時間軸方向には ∆t

u(x, t)

u(x + ∆x, y)

i j

u(x, t) = φ(i∆x, j∆t)

= u i j

(38)

と表現を改める。このようにすると、式

(37)

u i j+1 = 2u i j − u i j−1 + α (u i+1 j − 2u i j + u i−1 j ) (39)

α = ∆t ²

∆x ² (40)

である。