応用数値解析特論第 4 回

(1)

応用数値解析特論第 4 ^回

〜1^次元Poisson方程式に対する有限要素法〜

かつらだ

桂田祐史^{まさし}

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

2020

年

10

月

12

日

かつらだ桂田

まさし

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

(2)

と表すことができるが、u_i を並べたu^∗= (u1,u2,· · ·,um)^⊤ について、連立1次方程式Au=f が得られることは原理的に分かっている。しかし、A= (⟨ϕj, ϕi⟩) や f を実際に計算するのは、やり方を知らないと案外難しい。有限要素法ではこのあたりが良く整備されていて、明快なアルゴリズムが確立されている。

有限要素 ek ごとに、要素係数行列、要素自由項ベクトルというものを求め、それから Aと f を “組み立てる”。後半の操作を構造力学の用語にちなみ直接剛性法と呼ぶ。

(21)

4.3 有限要素解を求めるアルゴリズム

前項までに有限要素解uˆは定義された。uˆは ˆ

u=αϕ0+ Xm

i=1

uiϕi

と表すことができるが、u_i を並べたu^∗= (u1,u2,· · ·,um)^⊤ について、連立1次方程式Au=f が得られることは原理的に分かっている。しかし、A= (⟨ϕj, ϕi⟩) や f を実際に計算するのは、やり方を知らないと案外難しい。有限要素法ではこのあたりが良く整備されていて、明快なアルゴリズムが確立されている。

有限要素 ek ごとに、要素係数行列、要素自由項ベクトルというものを求め、それからAと f を “組み立てる”。後半の操作を構造力学の用語にちなみ直接剛性法と呼ぶ。

かつらだ桂田

まさし

(22)

4.3.1 長さ座標

各要素ek において

L₀(x_k₋₁) = 1, L₀(x_k) = 0, L1(xk−1) = 0, L1(xk) = 1 で定まる1次関数L₀,L₁を e_k の長さ座標と呼ぶ。

(6) L₀+L₁≡1

が成り立つ。また節点の座標xk−1,xk を用いて

(7) L0(x) = xk−x

x_k −x_k₋₁, L1(x) = x−xk−1

x_k−x_k₋₁ と具体的に表わせる (こちらはϕ_i と違って、後で用いる)。

ˆ

w ∈Xe に対してw_i := ˆw(x_i)とおくと、次式が成り立つ: (8) wˆ(x) =wk−1L0(x) +wkL1(x) =

X1 j=0

wk+j−1Lj (x ∈ek).

(23)

4.3.1 長さ座標

(6) L0+L1≡1

が成り立つ。

また節点の座標xk−1,xk を用いて

(7) L0(x) = xk−x

x_k −x_k₋₁, L1(x) = x−xk−1

ˆ

X1 j=0

かつらだ桂田

まさし

(24)

4.3.1 長さ座標

(6) L0+L1≡1

(7) L0(x) = xk−x

x_k −x_k₋₁, L1(x) = x−xk−1

ˆ

X1 j=0

(25)

4.3.1 長さ座標

(6) L0+L1≡1

(7) L0(x) = xk−x

x_k −x_k₋₁, L1(x) = x−xk−1

ˆ

w ∈Xe に対してw_i := ˆw(x_i)とおくと、次式が成り立つ:

(8) wˆ(x) =wk−1L0(x) +wkL1(x) = X1

j=0

かつらだ桂田

まさし

(26)

4.3.2 弱形式の分割

各要素ek について (9) ⟨u,v⟩_e_k :=

Z xk

xk−1

u^′(x)v^′(x)dx, (f,v)e_k :=

Z xk

xk−1

f(x)v(x)dx とおくと、

Galerkin法の弱形式

⟨u,ˆ vˆ⟩= (f,vˆ) +βvˆ(1) (ˆv∈Xˆ) は

(10)

Xm k=1

⟨u,ˆ vˆ⟩e_k = Xm k=1

(f,vˆ)e_k+βvˆ(1) (ˆv ∈Xˆ) と書き直せる。

(27)

4.3.2 弱形式の分割

各要素ek について (9) ⟨u,v⟩_e_k :=

Z xk

xk−1

u^′(x)v^′(x)dx, (f,v)e_k :=

Z xk

xk−1

f(x)v(x)dx とおくと、Galerkin法の弱形式

⟨u,ˆ vˆ⟩= (f,vˆ) +βvˆ(1) (ˆv∈Xˆ) は

(10)

Xm k=1

⟨u,ˆ vˆ⟩e_k = Xm k=1

(f,vˆ)e_k+βvˆ(1) (ˆv ∈Xˆ) と書き直せる。

かつらだ桂田

まさし

(28)

4.3.3 要素係数行列 , 要素自由項ベクトル

⟨u,ˆ vˆ⟩_e_k =

* ₁ X

j=0

u_k_+j₋₁L_j, X1

i=0

v_k_+i₋₁L_i +

e_k

= X1

j=0

X1 i=0

u_k+j₋₁v_k_+i₋₁⟨L_j,L_i⟩_e_k

= X1

i=0

X1 j=0

vk+i−1A^(k_ij⁾uk+j−1, ただし

A^(k)_ij :=⟨Lj,Li⟩_e_k. 一方、

(f,v)ˆ _e_k = (f, X1

j=0

v_k_+j₋₁L_j)_e_k = X1

j=0

v_k_+j₋₁(f,L_j)_e_k = X1

j=0

v_k_+j₋₁f_j^(k), ただし

f_j^(k):= (f,Lj)e_k.

(29)

4.3.3 要素係数行列 , 要素自由項ベクトル

そこで

uk :=

uk−1

uk

, vk :=

vk−1

vk

, fk := f₀^(k) f₁^(k)

! , (11a)

Ak := A^(k)₀₀ A^(k)₀₁ A^(k)₁₀ A^(k)₁₁

! , (11b)

gm:=

0 β

(11c)

とおくと、次式が得られる。

(12) ⟨u,ˆ vˆ⟩_e_k =v_k^⊤A_kuk, (f,v)ˆ _e_k =v_k^⊤fk (k = 1,· · ·,m), βv(1) =ˆ v_m^⊤gm. ここで^⊤ は転置(transpose)を表す。b^⊤a=a·b である。

uk,vk は要素節点パラメーター・ベクトル、fk は要素自由項ベクトル、A_k は要 素係数行列と呼ばれる。

かつらだ桂田

まさし

(30)

4.3.3 要素係数行列 , 要素自由項ベクトル

そこで

uk :=

uk−1

uk

, vk :=

vk−1

vk

, fk := f₀^(k) f₁^(k)

! , (11a)

Ak := A^(k)₀₀ A^(k)₀₁ A^(k)₁₀ A^(k)₁₁

! , (11b)

gm:=

0 β

(11c)

とおくと、次式が得られる。

(12) ⟨u,ˆ vˆ⟩_e_k =v_k^⊤A_kuk, (f,v)ˆ _e_k =v_k^⊤fk (k = 1,· · ·,m), βv(1) =ˆ v_m^⊤gm. ここで^⊤ は転置(transpose)を表す。b^⊤a=a·b である。

u ,v は要素節点パラメーター・ベクトル、f は要素自由項ベクトル、A は要

(31)

4.3.3 要素係数行列 , 要素自由項ベクトル

後のために実際にA_k, fk を求めよう。

⟨L_i,L_j⟩_e_k = Z x_k

x_k−1

L^′_i(x)L^′_j(x)dx= Z x_k

x_k−1

ε

(xk−xk−1)²dx= ε xk−xk−1

,

ε:=

(

1 (i=j)

−1 (i̸=j) であるから、

Ak = 1 x_k −x_k₋₁

1 −1

−1 1

.

一方

f_j^(k)= (f,Lj)e_k = Z xk

x_k−1

f(x)Lj(x)dx (j= 0,1).

この右辺の積分はf に応じて何らかの手段で計算しておく。

かつらだ桂田

まさし

(32)

4.3.3 要素係数行列 , 要素自由項ベクトル

(このスライドは最初はスキップしても良い。)

f が複雑な関数の場合は、f_j^(k) は厳密に計算できないかもしれないが、節点での値さえ分かれば近似計算は難しくない。

f ≒f(xk−1)L0+f(xk)L1 (onek) と1次補間近似を利用して

f_j^(k⁾= (f,Lj)e_k ≒(f(xk−1)L0+f(xk)L1,Lj)e_k =f(xk−1)(L0,Lj)e_k+f(xk)(L1,Lj)e_k. ところで(例えばx=x_k₋₁+ (x_k−x_k₋₁)t (0≤t ≤1)と変数変換して)

(L₀,L₀)_e_k = (L₁,L₁)_e_k = (x_k−x_k₋₁) Z 1

0

t²dt= x_k−x_k₋₁

3 ,

(L₀,L₁)_e_k = (L₀,L₁)_e_k = (x_k−x_k₋₁) Z 1

0

t(1−t)dt=x_k −x_k₋₁ 6

であるから

(33)

以下の話で必要になる式を再掲しておく。

弱形式は次のように書き直される。

(13)

Xm k=1

⟨u,

ˆ

v

ˆ

⟩_e_k

=

Xm k=1

(f

,v)

ˆ

_e_k

+

βv(1)

ˆ (ˆ

v ∈X

ˆ )

u_k

=

u_k₋₁

u_k

,

v_k

=

v_k₋₁

v_k

,

さらにf_k

,

A_k

,

g_m を適当に定義すると

⟨u,

ˆ

v

ˆ

⟩_e_k

=

v_k^⊤A_ku_k,

(f

,v)

ˆ

_e_k

=

v_k^⊤f_k

(k = 1,

· · · ,m), βv

ˆ (1) =

v_m^⊤g_m.

かつらだ桂田

まさし

(34)

4.3.4 直接剛性法 ( 近似方程式の組み立て )

(11a), (11b), (11c)で与えたベクトル、行列をm+ 1^{次元に拡大する。まず}

u:=



 u0

.. . um



, v :=



 v0

.. . vm





とおく。繰り返しになるが、ui = ˆu(xi),vi = ˆv(xi).

fk,Ak,g_m^∗については、0^{を補って、}R^m+1,M(m+ 1;R)^{の元に拡大する}:

f_k^∗:=





 0 .. . 0 f₀^(k) f₁^(k) 0 .. . 0







, A^∗_k:=







0 0 0

0

^A^(k)⁰⁰ ^A^(k)⁰¹

A^(k)₁₀ A^(k)₁₁

0 0 0 0





 (k= 1,· · ·,m), g_m^∗:=





 0

.. . 0 β





.

(35)

4.3.4 直接剛性法 ( 近似方程式の組み立て )

弱形式を書き直した(13)は次のように書き直される。

(15)

Xm k=1

v^⊤A^∗ku= Xm k=1

v^⊤f_k^∗+v^⊤g_m^∗ (v ∈Y).

ここでY は、vˆ= Xm

i=0

viϕi がXˆ に属するようなv = (v0,· · ·,vm)^⊤の全体、すなわち Y =

n

(v0,· · ·,vm)^⊤∈R^m+1v0= 0 o

. (15)は

v^⊤ Xm k=1

A^∗_k

! u=v^⊤

Xm k=1

f_k^∗+g_m^∗

!

(v ∈Y) と書き直せる。ゆえに

(16) A^∗:=

Xm k=1

A^∗_k, f^∗:= Xm k=1

f_k^∗+g_m^∗

とおけば

v^⊤A^∗u=v^⊤f^∗ (v ∈Y). すなわち

(17) v^⊤(A^∗u−f^∗) = 0 (v ∈Y).

かつらだ桂田

まさし

(36)

4.3.4 直接剛性法 ( 近似方程式の組み立て )

(15)

Xm k=1

v^⊤A^∗ku= Xm k=1

v^⊤f_k^∗+v^⊤g_m^∗ (v ∈Y).

i=0

n

(v0,· · ·,vm)^⊤∈R^m+1v0= 0 o

.

(15)は

v^⊤ Xm k=1

A^∗_k

! u=v^⊤

Xm k=1

f_k^∗+g_m^∗

!

(16) A^∗:=

Xm k=1

A^∗_k, f^∗:= Xm k=1

f_k^∗+g_m^∗

とおけば

(17) v^⊤(A^∗u−f^∗) = 0 (v ∈Y).

(37)

4.3.4 直接剛性法 ( 近似方程式の組み立て )

(15)

Xm k=1

v^⊤A^∗ku= Xm k=1

v^⊤f_k^∗+v^⊤g_m^∗ (v ∈Y).

i=0

n

(v0,· · ·,vm)^⊤∈R^m+1v0= 0 o

. (15)は

v^⊤ Xm k=1

A^∗_k

! u=v^⊤

Xm k=1

f_k^∗+g_m^∗

!

(v ∈Y) と書き直せる。

ゆえに

(16) A^∗:=

Xm k=1

A^∗_k, f^∗:= Xm k=1

f_k^∗+g_m^∗

とおけば

(17) v^⊤(A^∗u−f^∗) = 0 (v ∈Y).

かつらだ桂田

まさし

(38)

4.3.4 直接剛性法 ( 近似方程式の組み立て )

(15)

Xm k=1

v^⊤A^∗ku= Xm k=1

v^⊤f_k^∗+v^⊤g_m^∗ (v ∈Y).

i=0

n

(v0,· · ·,vm)^⊤∈R^m+1v0= 0 o

. (15)は

v^⊤ Xm k=1

A^∗_k

! u=v^⊤

Xm k=1

f_k^∗+g_m^∗

!

(16) A^∗:=

Xm k=1

A^∗_k, f^∗:=

Xm k=1

f_k^∗+g_m^∗ とおけば

v^⊤A^∗u=v^⊤f^∗ (v ∈Y).

(39)

4.3.4 直接剛性法 ( 近似方程式の組み立て )

(再掲17) v^⊤(A^∗u−f^∗) = 0 (v ∈Y).

これは次と同値である。

A^∗u−f ∈Y^⊥= n

(λ,0,· · ·,0)^⊤∈R^m+1λ∈Ro . つまりは

(A^∗u−f^∗)の最初の成分以外= 0. すなわち

(18) A^∗∗u=f^∗∗.

ここで

A^∗∗:=A^∗の第0行を除いたm×(m+ 1)行列, f^∗∗:=f^∗の第0成分を除いたm次元縦ベクトル.

(この段階で、係数行列A^∗∗ は正方行列ではない。しかしu の成分のうちu0は分かっていて未知ではない。その部分を右辺に移項することができ、そうすると正方行列係数の方程式が出来る)

部分配列を表すためのMATLAB風の記法を使うと、A^∗∗=A^∗(1 :m,0 :m), f^∗∗=f^∗(1 :m)と書ける。この記法は便利なので以下でも使うことにする。

かつらだ桂田

まさし

応用数値解析特論 第 4 回