応用数値解析特論第 6 回

(1)

応用数値解析特論第 6 回

〜

2

^次元

Poisson

方程式に対する有限要素法

(2)

かつらだ

桂田祐史^{まさし}

http://nalab.mind.meiji.ac.jp/~mk/ana/

2020

年

10

月

26

日

かつらだ桂田

まさし

祐史 http://nalab.mind.meiji.ac.jp/~mk/ana/応用数値解析特論第6回 2020年10月26日 1 / 24

(2)

本日の内容・今後の予定など

菊地

[1]

第

5

章に従い、空間

2

次元の有限要素法

(2

回目

)

で直接剛性法の解説を行なう。

菊地

[1]

第

7

章に掲載されていたサンプル

FORTRAN

プログラムを、

C

言語に翻訳したもの紹介する。実際に動作させて説明するのは、次回になる。

次回

(

来週は大学祭週間なので

11

^月

9

^日

)

^{からしばらく}

FreeFem++

を用いた数値実験を行なう。

(

スライドで説明しているため、例年よりも早く進行できると考えていましたが、実際にやってみると意外に時間のかかる面もあり、結局は例年とあまり変わらなくなっています。座学が長く続いてしまったことは、反省しています。

)

かつらだ桂田

まさし

(4)

本日の内容・今後の予定など

菊地

[1]

第

5

2 (2

回目

)

菊地

[1]

第

7 FORTRAN

C

次回

(

11

^月

9

^日

)

FreeFem++

(

)

かつらだ桂田

まさし

(5)

本日の内容・今後の予定など

菊地

[1]

第

5

2 (2

回目

)

菊地

[1]

第

7 FORTRAN

C

次回

(

11

^月

9

^日

)

FreeFem++

(

)

かつらだ桂田

まさし

(6)

本日の内容・今後の予定など

菊地

[1]

第

5

2 (2

回目

)

菊地

[1]

第

7 FORTRAN

C

次回

(

11

^月

9

^日

)

FreeFem++

(

)

かつらだ桂田

まさし

(7)

本日の内容・今後の予定など

菊地

[1]

第

5

2 (2

回目

)

菊地

[1]

第

7 FORTRAN

C

次回

(

11

^月

9

^日

)

FreeFem++

(

スライドで説明しているため、例年よりも早く進行できると考えていましたが、実際にやってみると意外に時間のかかる面もあり、結局は例年とあまり変わらなくなっています。座学が長く続いてしまったことは、

反省しています。

)

かつらだ桂田

まさし

(8)

前回にどこまで出来ていたか思い出す

弱形式は次の方程式と同値である。

(5.10)

N_e

X

k=1

⟨bu,bv⟩_e_k =

N_e

X

k=1

(f,v)be_k+ [g2,vb] (bv∈Xb).

本日の取り決め: (時間の都合上)g₂= 0とする。(5.10)の右辺第2項= 0.

任意の要素ek に注目し、e_k に含まれる節点に反時計回りにN0,N1,N2 と番号

(局所節点番号)をつけ

uⁱ = ˆu(Ni), vⁱ= ˆv(Ni) (i= 0,1,2),

uk = (u⁰u¹u²)^⊤, vk = (v⁰v¹v²)^⊤ とおいた。

3次元ベクトル fk, 3次正方行列Ak で

(5.13) ⟨u,ˆ vˆ⟩e_k =v_k^⊤A_kuk, (f,v)ˆ _e_k =v_k^⊤fk

を満たすものを求めた (要素自由項ベクトル,要素係数行列と呼ぶことにした)。

かつらだ桂田

まさし

(9)

前回にどこまで出来ていたか思い出す

(5.10)

N_e

X

k=1

⟨bu,bv⟩_e_k =

N_e

X

k=1

任意の要素ek に注目し、e_k に含まれる節点に反時計回りに N0,N1,N2 と番号

uⁱ = ˆu(Ni), vⁱ= ˆv(Ni) (i= 0,1,2),

uk = (u⁰u¹u²)^⊤, vk = (v⁰v¹v²)^⊤ とおいた。

(5.13) ⟨u,ˆ vˆ⟩e_k =v_k^⊤A_kuk, (f,v)ˆ _e_k =v_k^⊤fk

かつらだ桂田

まさし

(10)

前回にどこまで出来ていたか思い出す

(5.10)

N_e

X

k=1

⟨bu,bv⟩_e_k =

N_e

X

k=1

任意の要素ek に注目し、e_k に含まれる節点に反時計回りに N0,N1,N2 と番号

uⁱ = ˆu(Ni), vⁱ= ˆv(Ni) (i= 0,1,2),

uk = (u⁰u¹u²)^⊤, vk = (v⁰v¹v²)^⊤ とおいた。

(5.13) ⟨u,ˆ vˆ⟩e_k =v_k^⊤Akuk, (f,v)ˆ e_k =v_k^⊤fk

かつらだ桂田

まさし

(11)

5.4 近似方程式の組み立て — 直接剛性法行列・ベクトルの拡大

前小節(§5.3)で与えたベクトル、行列をm次元に拡大する。まず、

ui :=u(Pb i),vi :=bv(Pi) (i= 1,2,· · · ,m)として、次のようにおく。

u :=





 u1

u2

... u_m





, v :=





 v1

v2

... v_m





.

有限要素ek (k = 1,2,· · ·,Ne)の節点N0,N1,N2に対して、 N₀=P_i_k,0, N₁=P_i_k,1, N₂=P_i_k,2

となる整数ik,0,ik,2,ik,2 を取る(これらを全体節点番号と呼ぶ)。このときuk の成分u⁰=u(Nb ₀),u¹=u(Nb ₁),u²=u(Nb ₂)について

u⁰=ui_k,0, u¹=ui_k,1, u²=ui_k,2

が成り立つ。

かつらだ桂田

まさし

(12)

5.4 近似方程式の組み立て — 直接剛性法行列・ベクトルの拡大

前小節(§5.3)で与えたベクトル、行列をm次元に拡大する。まず、

ui :=u(Pb i),vi :=bv(Pi) (i= 1,2,· · · ,m)として、次のようにおく。

u :=





 u1

u2

... u_m





, v :=





 v1

v2

... v_m





.

有限要素ek (k = 1,2,· · ·,Ne)の節点N0,N1,N2に対して、

N₀=P_i_k,0, N₁=P_i_k,1, N₂=P_i_k,2

となる整数ik,0,ik,2,ik,2 を取る(これらを全体節点番号と呼ぶ)。このときuk の成分u⁰=u(Nb ₀),u¹=u(Nb ₁),u²=u(Nb ₂)について

u⁰=ui_k,0, u¹=ui_k,1, u²=ui_k,2

が成り立つ。

かつらだ桂田

まさし

(13)

5.4 直接剛性法行列・ベクトルの拡大

f

_k

= (f

₀^(k⁾

f

₁^(k)

f

₂^(k)

)

^⊤ ^{である。これを}

m

^{次元ベクトル}

f

_k^∗ に拡大する。

f

_k^∗

∈ R

^m ^を

i

k,0 成分

= f

₀^(k)

, i

k,1 成分

= f

₁^(k)

, i

k,2 成分

= f

₂^(k) ^{で、それ以} 外の成分はすべて

0

であるようなベクトルとする。

例えば

i

₀^(k⁾

< i

₁^(k)

< i

₂^(k) ならば

1 i

_k_,0

i

_k,1

i

_k,2

m

↓ ↓ ↓ ↓ ↓

f

_k^∗

= (0 · · · 0 f

₀^(k)

0 · · · 0 f

₁^(k)

0 · · · 0 f

₂^(k)

0 · · · 0)

^⊤

.

このように定義すると、

(1) (f , b v)

ek

= v

_k^⊤

f

_k

= v

^⊤

f

_k^∗

(k = 1, 2, · · · , N

e

).

かつらだ桂田

まさし

(14)

5.4 直接剛性法行列・ベクトルの拡大

f

_k

= (f

₀^(k⁾

f

₁^(k)

f

₂^(k)

)

^⊤ ^{である。これを}

m

^{次元ベクトル}

f

_k^∗ に拡大する。

f

_k^∗

∈ R

^m ^を

i

k,0 成分

= f

₀^(k)

, i

k,1 成分

= f

₁^(k)

, i

k,2 成分

= f

₂^(k) ^{で、それ以} 外の成分はすべて

0

であるようなベクトルとする。

例えば

i

₀^(k⁾

< i

₁^(k)

< i

₂^(k) ならば

1 i

_k_,0

i

_k,1

i

_k,2

m

↓ ↓ ↓ ↓ ↓

f

_k^∗

= (0 · · · 0 f

₀^(k)

0 · · · 0 f

₁^(k)

0 · · · 0 f

₂^(k)

0 · · · 0)

^⊤

.

このように定義すると、

(1) (f , b v)

ek

= v

_k^⊤

f

_k

= v

^⊤

f

_k^∗

(k = 1, 2, · · · , N

e

).

かつらだ桂田

まさし

(15)

5.4 直接剛性法行列・ベクトルの拡大

同様の考え方で、3^{次正方行列}Ak= (A^(k)_ij )^を、m^{次正方行列}A^∗k= (a^(k)_ij )^に a_i^(k)

k,0i_k,0=A^(k)₀₀, a^(k)_i

k,0i_k,1=A^(k)₀₁, a^(k)_i

k,0i_k,2=A^(k)₀₂, a_i^(k)

k,1i_k,0=A^(k)₁₀, a^(k)_i

k,1i_k,1=A^(k)₁₁, a^(k)_i

k,1i_k,2=A^(k)₁₂, a_i^(k)

k,2i_k,0=A^(k)₂₀, a^(k)_i

k,2i_k,1=A^(k)₂₁, a^(k)_i

k,2i_k,2=A^(k)₂₂, それ以外= 0

として拡大する。

例えばik,0<ik,1<ik,2ならば

A^∗k=







ik,0 ik,1 ik,2

↓ ↓ ↓

A^(k)₀₀ A^(k)₀₁ A^(k)₀₂ ←ik,0

A^(k)₁₀ A^(k)₁₁ A^(k)₁₂ ←ik,1

A^(k)₂₀ A^(k)₂₁ A^(k)₂₂ ←ik,2







(書いてないとこは0).

これを用いると

(2) ⟨u,bvb⟩e_k =v_k^⊤Akuk=v^⊤A^∗_ku (k= 1,2,· · ·,Ne).

かつらだ桂田

まさし

(16)

5.4 直接剛性法行列・ベクトルの拡大

k,0i_k,0=A^(k)₀₀, a^(k)_i

k,0i_k,1=A^(k)₀₁, a^(k)_i

k,0i_k,2=A^(k)₀₂, a_i^(k)

k,1i_k,0=A^(k)₁₀, a^(k)_i

k,1i_k,1=A^(k)₁₁, a^(k)_i

k,1i_k,2=A^(k)₁₂, a_i^(k)

k,2i_k,0=A^(k)₂₀, a^(k)_i

k,2i_k,1=A^(k)₂₁, a^(k)_i

k,2i_k,2=A^(k)₂₂, それ以外= 0

として拡大する。例えばik,0<ik,1<ik,2ならば

A^∗k=







ik,0 ik,1 ik,2

↓ ↓ ↓

A^(k)₀₀ A^(k)₀₁ A^(k)₀₂ ←ik,0

A^(k)₁₀ A^(k)₁₁ A^(k)₁₂ ←ik,1

A^(k)₂₀ A^(k)₂₁ A^(k)₂₂ ←ik,2







かつらだ桂田

まさし

(17)

5.4 直接剛性法行列・ベクトルの拡大

k,0i_k,0=A^(k)₀₀, a^(k)_i

k,0i_k,1=A^(k)₀₁, a^(k)_i

k,0i_k,2=A^(k)₀₂, a_i^(k)

k,1i_k,0=A^(k)₁₀, a^(k)_i

k,1i_k,1=A^(k)₁₁, a^(k)_i

k,1i_k,2=A^(k)₁₂, a_i^(k)

k,2i_k,0=A^(k)₂₀, a^(k)_i

k,2i_k,1=A^(k)₂₁, a^(k)_i

k,2i_k,2=A^(k)₂₂, それ以外= 0

として拡大する。例えばik,0<ik,1<ik,2ならば

A^∗k=







ik,0 ik,1 ik,2

↓ ↓ ↓

A^(k)₀₀ A^(k)₀₁ A^(k)₀₂ ←ik,0

A^(k)₁₀ A^(k)₁₁ A^(k)₁₂ ←ik,1

A^(k)₂₀ A^(k)₂₁ A^(k)₂₂ ←ik,2







かつらだ桂田

まさし

(18)

5.4 直接剛性法全体係数行列と全体自由項ベクトル

ゆえに弱形式

(5.10)

は

Ne

X

k=1

v

^⊤

A

^∗_k

u =

Ne

X

k=1

v

^⊤

f

_k^∗

すなわち

v

^⊤

Ne

X

k=1

A

^∗_k

u −

Ne

X

k=1

f

_k^∗

!

= 0

と同値になる。

そこで

A

^∗

:=

Ne

X

k=1

A

^∗_k

, f

^∗

:=

Ne

X

k=1

f

_k^∗

とおけば

(3) v

^⊤

(A

^∗

u − f

^∗

) = 0.

かつらだ桂田

まさし

(19)

5.4 直接剛性法全体係数行列と全体自由項ベクトル

ゆえに弱形式

(5.10)

は

Ne

X

k=1

v

^⊤

A

^∗_k

u =

Ne

X

k=1

v

^⊤

f

_k^∗

すなわち

v

^⊤

Ne

X

k=1

A

^∗_k

u −

Ne

X

k=1

f

_k^∗

!

= 0

そこで

A

^∗

:=

Ne

X

k=1

A

^∗_k

, f

^∗

:=

Ne

X

k=1

f

_k^∗

とおけば

(3) v

^⊤

(A

^∗

u − f

^∗

) = 0.

かつらだ桂田

まさし

(20)

5.4 直接剛性法全体係数行列と全体自由項ベクトル

ゆえに弱形式

(5.10)

は

Ne

X

k=1

v

^⊤

A

^∗_k

u =

Ne

X

k=1

v

^⊤

f

_k^∗

すなわち

v

^⊤

Ne

X

k=1

A

^∗_k

u −

Ne

X

k=1

f

_k^∗

!

= 0

そこで

A

^∗

:=

Ne

X

k=1

A

^∗_k

, f

^∗

:=

Ne

X

k=1

f

_k^∗

とおけば

(3) v

^⊤

(A

^∗

u − f

^∗

) = 0.

かつらだ桂田

まさし

(21)

5.4 直接剛性法全体係数行列と全体自由項ベクトル

ここで

v

は

Y := {



  v

1

.. . v

m



  ∈ R

^m

; v

_i

= 0 (P

_i

∈ Γ

₁ なる

i) }

の任意の元であるから、

(3)

は次と同値である。

A

^∗

u − f ∈ Y

^⊥

= {



  w

1

.. . w

_m



  ∈ R

^m

; w

_i

= 0 (P

_i

̸∈ Γ

₁ なる

i) }

すなわち

(4) A

^∗∗

u = f

^∗∗

.

ここで

A

^∗∗

:= A

^∗ の第

i

^行

(P

i

∈ Γ

1 なる

i)

^{を除いた行列}

,

f

^∗∗

:= f

^∗ の第

i

成分

(P

_i

∈ Γ

₁ なる

i)

を除いた縦ベクトル

.

ところで

u

_i

= g

₁

(P

_i

) (P

_i

∈ Γ

₁ なる

i )

があるから、これを代入して消去できる。

A

^∗∗ ^{を列ベクトルで}

A

^∗∗

= (a

₁

a

₂

· · · a

_m

)

のように表示すると、

(4)

^は

X

m

i=1

a

_i

u

_i

= f

^∗∗

.

この左辺を

X

Pi∈Γ1なるi

a

_i

u

i

+ X

Pi̸∈Γ1なるi

a

_i

u

i

= f

^∗∗ と分解して、移項すると

X

Pi ̸∈Γ1なるi

a

_i

u

i

= f

^∗∗

− X

Pi∈Γ1なるi

a

_i

u

i

.

かつらだ桂田

まさし

(22)

5.4 直接剛性法全体係数行列と全体自由項ベクトル

ここで

v

は

Y := {



  v

1

.. . v

m



  ∈ R

^m

; v

_i

= 0 (P

_i

∈ Γ

₁ なる

i) }

(3)

A

^∗

u − f ∈ Y

^⊥

= {



  w

1

.. . w

_m



  ∈ R

^m

; w

_i

= 0 (P

_i

̸∈ Γ

₁ なる

i) }

すなわち

(4) A

^∗∗

u = f

^∗∗

.

ここで

A

^∗∗

:= A

^∗ の第

i

^行

(P

i

∈ Γ

1 なる

i)

,

f

^∗∗

:= f

^∗ の第

i

成分

(P

_i

∈ Γ

₁ なる

i)

.

ところで

u

_i

= g

₁

(P

_i

) (P

_i

∈ Γ

₁ なる

i )

A

^∗∗ ^{を列ベクトルで}

A

^∗∗

= (a

₁

a

₂

· · · a

_m

)

(4)

^は

X

m

i=1

a

_i

u

_i

= f

^∗∗

.

この左辺を

X

Pi∈Γ1なるi

a

_i

u

i

+ X

Pi̸∈Γ1なるi

a

_i

u

i

= f

^∗∗ と分解して、移項すると

X

Pi ̸∈Γ1なるi

a

_i

u

i

= f

^∗∗

− X

Pi∈Γ1なるi

a

_i

u

i

.

かつらだ桂田

まさし

(23)

5.4 直接剛性法全体係数行列と全体自由項ベクトル

ここで

v

は

Y := {



  v

1

.. . v

m



  ∈ R

^m

; v

_i

= 0 (P

_i

∈ Γ

₁ なる

i) }

(3)

A

^∗

u − f ∈ Y

^⊥

= {



  w

1

.. . w

_m



  ∈ R

^m

; w

_i

= 0 (P

_i

̸∈ Γ

₁ なる

i) }

すなわち

(4) A

^∗∗

u = f

^∗∗

.

ここで

A

^∗∗

:= A

^∗ の第

i

^行

(P

i

∈ Γ

1 なる

i )

,

f

^∗∗

:= f

^∗ の第

i

成分

(P

_i

∈ Γ

₁ なる

i)

.

ところで

u

_i

= g

₁

(P

_i

) (P

_i

∈ Γ

₁ なる

i )

A

^∗∗ を列ベクトルで

A

^∗∗

= (a

₁

a

₂

· · · a

_m

)

(4)

^は

X

m

i=1

a

_i

u

_i

= f

^∗∗

.

この左辺を

X

Pi∈Γ1なるi

a

_i

u

i

+ X

Pi ̸∈Γ1なるi

a

_i

u

i

= f

^∗∗

と分解して、移項すると

X

Pi ̸∈Γ1なるi

a

_i

u

i

= f

^∗∗

− X

Pi∈Γ1なるi

a

_i

u

i

.

かつらだ桂田

まさし

(24)

5.4 直接剛性法全体係数行列と全体自由項ベクトル

ゆえに

Au

^∗

= f ,

ただし

u

^∗

:= u

の第

i

^成分

(P

i

∈ Γ

1 なる

i )

^{を除いた縦ベクトル}

, A := A

^∗∗^の第

i

^列

(P

i

∈ Γ

1 なる

i )

^{を除いた正方行列}

, f := f

^∗∗

− X

Pi ∈Γ1なるi

a

_i

u

_i

.

(P

_i

∈ Γ

₁ となる

i

について

u

_i

= ˆ u(P

_i

) = g

₁

(P

_i

)

は既知であり、未知数に含める必要はない。

)

かつらだ桂田

まさし

(25)

5.5 連立 1 次方程式の具体例

簡単な問題に対する有限要素法の連立

1

次方程式を実際に求めてみよう。

Ω = (0, 1) × (0, 1),

Γ

₁

= { (x, y); x = 0, 0 ≤ y ≤ 1 } ∪ { (x, y); 0 ≤ x ≤ 1, y = 0 } , Γ

₂

= { (x, y); x = 1, 0 < y ≤ 1 } ∪ { (x, y); 0 < x ≤ 1, y = 1 } , g

1

≡ 0, g

2

≡ 0, f ≡

^定数関数

f .

すなわち

−△ u = f (in Ω), u = 0 on Γ

₁

, ∂u

∂n = 0 on Γ

₂

.

かつらだ桂田

まさし

(26)

5.5 連立 1 次方程式の具体例

簡単な問題に対する有限要素法の連立

1

次方程式を実際に求めてみよう。

Ω = (0, 1) × (0, 1),

Γ

₁

= { (x, y); x = 0, 0 ≤ y ≤ 1 } ∪ { (x, y); 0 ≤ x ≤ 1, y = 0 } , Γ

₂

= { (x, y); x = 1, 0 < y ≤ 1 } ∪ { (x, y); 0 < x ≤ 1, y = 1 } , g

1

≡ 0, g

2

≡ 0, f ≡

^定数関数

f .

すなわち

−△ u = f (in Ω), u = 0 on Γ

₁

, ∂u

∂n = 0 on Γ

₂

.

かつらだ桂田

まさし

(27)

5.5 連立 1 次方程式の具体例

- 6

O x

y

x = 1 y = 1

Ω

図1:領域Ω

- 6

x y

P

0

P

1

P

2

P

3

P

₄

P

₅

P

₆

P

7

P

8

e

₀

e

₁

e

₂

e

₃

e

₄

e

₅

e

₆

e

₇

図 2:要素分割

かつらだ桂田

まさし

(28)

5.5 連立 1 次方程式の具体例

正方形領域

Ω

を図

2

のように

3

角形要素によって要素分割する。

有限要素は次の二つのタイプがある

(

^タイプ

I, II

^{と呼ぶことにする}

)

^。各々に図

3

のように局所節点番号をつける

(

左下から反時計回り

)

。

タイプ

I e

0

, e

2

, e

4

, e

6

.

タイプ

II e

₁

, e

₃

, e

₅

, e

₇

.

N

0

N

1

N

2

Type I

N

0

N

1

N

₂

Type II

図3:二つのタイプの有限要素と局所節点番号

かつらだ桂田

まさし

(29)

5.5 連立 1 次方程式の具体例

タイプが同じならば、要素係数行列Ak,要素自由項ベクトルfk が等しいことはすぐ分かる。それぞれAI,AII,fI,fII で表すことにする。

タイプIについては

D =h², S = h² 2 ,

AI = 1 2



 1 −1 0

−1 2 −1

0 −1 1



, fI = f h² 6



 1 1 1



.

タイプIIについては

D =h², S = h² 2 ,

A_II = 1 2



 1 0 −1

0 1 −1

−1 −1 2



, fII =f h² 6



 1 1 1



.

かつらだ桂田

まさし

(30)

5.5 連立 1 次方程式の具体例

これらから、全体的な近似方程式を作ろう。

そのために局所的な節点番号と、全体的な節点番号の対応づけが必要である。そこで以下のような対応表を用意する。

要素 e0 e1 e2 e3 e4 e5 e6 e7

要素タイプ I II I II I II I II N0の全体節点番号 0 0 1 1 3 3 4 4 N1の全体節点番号 3 4 4 5 6 7 7 8 N2の全体節点番号 4 1 5 2 7 4 8 5

N0 N1

N2

Type I N0

N1

N2

Type II

- 6

x y

P0

P1

P2

P3

P4

P5

P6

P7

P8

e0

e1

e2

e3

e4

e5

e6

e7

かつらだ桂田

まさし

応用数値解析特論 第 6 回