応用数値解析特論第 6

(1)

応用数値解析特論第 6 回

〜

2

^次元

Poisson

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

(1)

かつらだ

桂田

まさし

祐史

https://m-katsurada.sakura.ne.jp/ana2022/

2023

^年

5

^月

23

^日

かつらだまさし

(2)

1

本日の内容など

2 1

次元

Poisson

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

(

続き

)

サンプル・プログラム

fem1d.c

問題

プログラムの解説実験

参考

:

昔の練習問題

3 2

次元の有限要素法

三角形要素への分割と区分的

1

次多項式三角形

e

上の

1

次関数

L

iと

(L

j

, L

i

)

e

,⟨L

j

, L

i

⟩

_e

三角形の面積

三角形要素の面積座標

L

i

面積座標の積の積分と

(L

j

, L

i

)

e

L

i

(x , y ) = a

i

+ b

i

x + c

i

y

の係数決定と

⟨ L

j

, L

i

⟩

_e 要素係数行列の計算

近似方程式の組み立て

—

直接剛性法連立

1

次方程式の具体例

プログラム

方程式を立てるのに必要なものサンプルプログラム紹介

4 C

言語による

2

次元有限要素法サンプル・プログラムの紹介進行表

試しに実行

有限要素解を求めるプログラム

naive, band

の理解プログラム

naive

の内部構造

参考課題

5 FreeFem++

の文法

はじめに

汎用のプログラミング機能

C

言語と良く似ているところデータ型

配列型

6

参考文献

(3)

本日の内容など

前回説明した

1

^次元

Poisson

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

C

^プログラムを読んでみる。

2

^次元の

Poisson

方程式に対する有限要素法で、連立１次方程式が具

体的にどのようになるか説明する。

(4)

4.5 サンプル・プログラム fem1d.c 4.5.1 問題

以下に紹介する

C

プログラム

fem1d.c

は

https://m-katsurada.sakura.ne.jp/program/fem/fem1d.c

に置いてある。現象数理学科

Mac

ならば、ターミナルから

curl -O https://m-katsurada.sakura.ne.jp/program/fem/fem1d.c

で入手できる。コンパイル、実行の仕方はプログラムの先頭部分に注釈として書いてある

(head fem1d.c

で先頭

10

行を見れば分かる)。

このプログラムが対象としている問題は、

f ≡ 1

で、境界条件は同次、すなわち

α = β = 0

の場合である。具体的に書き下すと

(1) − u

^′′

= 1, u(0) = u

^′

(1) = 0.

この問題の厳密解は

u(x) = x(2 − x)/2

である。

かつらだ桂田

まさし

祐史 https://m-katsurada.sakura.ne.jp/ana2022/応用数値解析特論第6回〜2次元Poisson方程式に対する有限要素法(1) 3 / 66

(5)

4.5.2 プログラムの解説

main()

を読むと分かるように、最初に

nnode

総節点数

(the number of nodes) nelmt

総要素数

(the number of elements)

nbc

ディリクレ境界にある接点の個数

(1

または

2) x[]

節点の座標

ibc

ディリクレ境界にある接点の節点番号を決めている。

連立

1

次方程式の構成は、関数

assem()

で行っている

(assemblage)

。その作業内容は

3

つに分かれる。

1

am, fm

を

0

クリアする。

2 すべての有限要素について、要素係数行列

ae,

要素自由ベクトル

fe

を関数

ecm()

で計算して

(element coeﬃcient matrix)、それぞれ全体

係数行列

am

、全体自由項ベクトル

fm

に算入する。

3 ディリクレ境界上にある節点に対応する部分を修正する。

(6)

4.5.2 プログラムの解説

関数

ecm()

で必要となる事項の復習。

e

k

= [x

k−1

, x

k

]

^{とすると、}

A

k

= 1 x

k

− x

k−1

1 −1

− 1 1

, f

k

=

(f , L

0

)

e_k

(f , L

1

)

e_k

であった。

f

を

1

次近似することにすれば、

f

_k

≒ x

k

− x

k−1

6

2f (x

k−1

) + f (x

k

) f (x

_k−1

) + 2f (x

k

)

.

前回の授業の説明では、

A

k を

m + 1

次正方行列

A

^∗_k に

“

拡大して

”

、

A

^∗

:=

X

m k=1

A

^∗k を計算し、その最初の列と最初の行

(Dirichlet

境界条件を課す

x

0に対応

)

を削除する。

ベクトル

f

k についても同様

(f

_k^∗

, f

^∗

:=

X

m k=1

f

_k^∗

+ g

_m^∗

)

。右辺に

α



  A

^∗₁₀

. . . A

^∗m0



 

^を移

項する

(α

は

Dirichlet

条件

u(0) = α

の

α)

。

この問題では

α = 0

なので移項せず第

0

列を

0

クリアするだけ。

列や行の削除は大変

(大規模コピー？)

なので、書き換えることに

かつらだする。

桂田まさし

(7)

4.5.3 実験

fem1d.c

のコンパイル＆実行＆

gnuplot

によるグラフ描画

$ cc -o fem1d fem1d.c

$ ./fem1d

nodal values of u (節点での u

の値)

i u i u i u

0 0.000e+00 1 9.500e-02 2 1.800e-01 3 2.550e-01 4 3.200e-01 5 3.750e-01 6 4.200e-01 7 4.550e-01 8 4.800e-01 9 4.950e-01 10 5.000e-01

$ cat fem1d.out 0.000000 0.000000 0.100000 0.095000 0.200000 0.180000 0.300000 0.255000 0.400000 0.320000 0.500000 0.375000 0.600000 0.420000 0.700000 0.455000 0.800000 0.480000 0.900000 0.495000 1.000000 0.500000

$ gnuplot

gnuplot> plot "fem1d.out" with lp, x*(2-x)/2

かつらだ桂田

まさし

(8)

4.5.3 実験

0 0.1 0.2 0.3 0.4 0.5

0 0.2 0.4 0.6 0.8 1

"fem1d.out"

x*(2-x)/2

図

1: fem1d.c

の計算結果

(m=10)

と厳密解 ^x(2⁻^x)

2 のグラフを重ね書き

かつらだ桂田

まさし

(9)

4.5.4 参考 : 昔の練習問題

FreeFem++

がまだなかった頃、有限要素法のプログラムを、

C

言語や

Fortran

のようなプログラミング言語でプログラムを書いていました。

そのときは

(

アルゴリズムの理解する助けになると考えて

)

以下のような練習問題を出していました。参考まで。

1 両側ディリクレ条件

u(0) = u(1) = 0

の問題を解く。

2 非同次ディリクレ条件

u(0) = α

3 非同次

Neumann

条件

u

^′

(1) = β

4

− (pu

^′

)

^′

= f

という一般の楕円型方程式の問題を解く。

(p

は

min

x

p(x) > 0

を満たす既知の関数)

(10)

5 2 次元の有限要素法

第

2

回の授業で扱った

Poisson

方程式の境界値問題を題材に、

2

次元領域における有限要素法を説明する

(

菊地

[?]

第

5

章

)

。

(

思い出す

) Ω

は

R

²の有界領域で、その境界

Γ

は区分的に十分滑らかであるとする。

また

Γ

1

, Γ

2は条件

Γ = Γ

1

∪ Γ

2

, Γ

1

∩ Γ

2

= ∅, Γ

1

̸= ∅

を満たすとする。

f : Ω → R , g

1

: Γ

1

→ R , g

2

: Γ

1

→ R

が与えられたとき、次の

Poisson

方程式の境界値問題を考える。

問題

(P)

次式を満たす

u

を求めよ

:

−△ u = f in Ω, (2)

u = g

1

on Γ

1

, (3)

∂ u

∂n = g

2

on Γ

2

, (4)

ここで

n

は

Γ

の外向き単位法線ベクトルを表す。

大筋は

1

次元の場合と同様だが、

(i)

各要素内の計算に面積座標を使うところと、

(ii)

直接剛性法が

2

次元でも実現可能であることを理解するところが要点となる。

かつらだ桂田

まさし

(11)

5.1 三角形要素への分割と区分的 1 次多項式

領域

Ω ⊂ R

²^に対し、

Ω

^を三角形

e

k

(1 ≤ k ≤ N

e

)

^{に分割する}

: Ω ≒ Ω := b

N_e

[

k=1

e

k

.

ただし、重なりや、すき間、頂点が他の三角形の辺上にあることは避けることにする。

[0]

[1]

[2]

[3]

[4]

[6][5]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17][18][19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[30][29]

[31]

[32]

[33]

[34]

[35]

[36]

[37]

[38]

[39]

[40]

[41][42][43]

[44]

[45]

[46]

[47]

[48]

[49]

[50]

[51]

[52]

[54] [53]

[55]

[56]

[57]

[58]

[59]

[60]

[61]

[62]

[63]

[64]

[65] [66] [67]

[68]

[69]

[70]

[71]

[72]

[73]

[74]

[75]

[76]

[78] [77]

[79]

[80]

[81]

[82]

[83]

[84]

[85]

[86]

[87]

[88]

[89] [90]

[91]

[92]

[93]

[94]

[95]

[96]

[97]

[98]

[99]

[100]

[101]

[102]

[103]

[104]

[105]

[106]

[107]

[108]

[109]

[110]

[111]

[112]

[113] [114]

[115]

[116]

[117]

[118]

[119]

[120]

[121]

[122]

[123]

[124]

[125]

[126]

[127]

[128]

[129]

[130]

[131]

[132]

[133]

[134]

[135]

[136]

[137] [138]

[139]

[140]

[141]

[142]

[143]

[144]

[145]

[146]

[147]

[148]

[149]

[150]

[151]

[152]

[153]

[154]

[155]

[156]

[157]

[158]

[159]

[160]

[161] [162]

[163]

[164]

[165]

[166]

[167]

[168]

[169]

[170]

[171]

[172]

[173]

[174]

[175]

[176]

[177]

[178]

[179]

[180]

[181]

[182]

[183]

[184]

[185] [186]

[187]

[188]

[189]

[190]

[191]

図

2:

二つの円で囲まれた閉領域

Ω

を三角形の合併で近似する

Ω

が多角形でない限り、境界は「曲がっている」。

Ω =

Ne

[

k=1

e

k は期待できない。

各三角形を三角形

(

有限

)

要素とよぶ。

(12)

5.1 三角形要素への分割と区分的 1 次多項式

N

e は要素の総数

(the number of elements)

で、プログラムでは

NELMT

のような名前の変数で記憶されることが多い。

有限要素の頂点を節点

(node)

と呼び、

{ P

_i

}

^mi=1 のように番号をつけておく。

m

は節点の総数

(the number of nodes)

で、プログラムでは

NNODE

のような名前の変数で記憶されることが多い。

注意

1

次元の場合、節点の個数

=

要素の個数

+ 1 (NELMT=m+1)

という簡単な関係が成立していた。(区間を

m

等分したとき、m を用いて、要素を

e

k

(1 ≤ k ≤ m),

節点を

x

k

(0 ≤ k ≤ m)

と番号づけることが出来た。)

2

次元の場合は、そのような関係はない。

かつらだ桂田

まさし

(13)

5.1 三角形要素への分割と区分的 1 次多項式

Ω b

上連続で、各有限要素

e

k 上で

x

と

y

の

1

次多項式関数に等しいものを、

区分的

1

次多項式と呼び、その全体を

X e

で表わす

( X e

はここだけの記号)。

{ e

k

}

^N_k₌₁^e と

X e

の対を区分的

1

次有限要素

(piecewise linear ﬁnite element)

と呼ぶ。

2

変数の

1

次関数

z = a + bx + cy

のグラフは平面であるから、

X e

のグラフは、空間内の三角形を連続につなげた

“折れ面”

である。

試行関数の空間

X ˆ

g₁

,

試験関数の空間

X ˆ

は次のように選ぶ。

(5) X ˆ

g₁

= n

ˆ

w ∈ X e w ˆ = ˆ g

1

on Γ

1

o

, X ˆ = n

ˆ

v ∈ X e v ˆ = 0 on Γ

1

o . (ˆ g

₁ は

g

₁に近い計算しやすい関数)

X e

の基底関数

{ ϕ

_i

}

^m_i=1は

ϕ

_i

∈ X e (i = 1, 2, · · · , m), (6a)

ϕ

i

(P

j

) = δ

ij

(i, j = 1, 2, · · · , m) (6b)

満たすものを採用する

(この条件で一意に確定し、線形独立であることに注意)。

(14)

5.2 三角形 e 上の 1 次関数 L i と (L j , L i ) e , ⟨ L j , L i ⟩ e

5.2.1

三角形の面積

平面上の三角形要素

e

を考える

(本来は、以下の N

i

, L

i も含めて、要素によるので、e_k

, N

^k_i

, L

^k_i のように書いた方が良いかもしれないが、うっとうしいので

k

は略する)。

e

に属する節点を、反時計まわりに

N

i

= (x

i

, y

i

) (i = 0, 1, 2)

とする。

後でしばしば必要になるので、e の面積

| e |

を計算しておこう。

| e | = 1 2 det

x

1

− x

0

x

2

− x

0

y

1

− y

0

y

2

− y

0

= 1

2 [(x

₁

− x

₀

)(y

₂

− y

₀

) − (y

₁

− y

₀

)(x

₂

− x

₀

)] . 1

次多項式全体を

P

¹ と表す:

P

¹

:= { u ˆ | ( ∃ α

₀

, α

₁

, α

₂

∈ R ) ˆ u(x, y) = α

₀

+ α

₁

x + α

₂

y } .

任意の

u ˆ ∈ P

¹は、3節点

N

_i における値

u

ⁱ

:= u(N b

_i

) (i = 0, 1, 2)

を指定すれば定まる

(u

_i

= u(P b

_i

)

と混同しないように上に添字をつけた)。これは、直観的に明らかであるが

(平面は、その上にある 3

点

(ただし同一直線上にはないとする)

を指定すれば定まる)、すぐ後で証明する。

かつらだ桂田

まさし

(15)

5.2.2 三角形要素の面積座標 L i

節点

N

i で

1,

他の節点で

0

となる

1

次関数を

L

i とする

(i = 0, 1, 2)。つまり L

_i

∈ P

¹

,

(7a)

L

i

(N

j

) = L

i

(x

j

, y

j

) = δ

ij

(j = 0, 1, 2).

(7b)

(L

0

, L

1

, L

2

)

は

P

¹の基底になる。実際、1次独立性はすぐ分かり、dim

P

¹

= 3

であるから、任意の

u b ∈ P

¹は、

(8) u(x, b y) =

X

2 i=0

u

ⁱ

L

_i

(x, y) ((x, y) ∈ e)

と表される。

任意の

P ∈ e

に対して、3 実数

(L

₀

(P), L

₁

(P), L

₂

(P))

を

P

の面積座標

(area coordinate)

あるいは重心座標

(barycentric coordinate)

と呼ぶ。

(色々裏がある

けれど今回は駆け足で進む。)

任意の

P ∈ e

に対して次式が成り立つ。

L

0

(P) + L

1

(P) + L

2

(P) = 1.

(P

¹の基底としては

3

つ必要だが、座標としては

2

つで十分ということになる。) 以下、L_i のグラフの鳥瞰図と等高線を描こう。

(16)

Mathematica のコード例と実行結果

xs = {0, 3, 1}; ys = {0, 1, 2};

{a, b, c} = Inverse[Transpose[{{1, 1, 1}, xs, ys}]];

L[i_, x_, y_] := a[[i]] + b[[i]]*x + c[[i]]*y

xmin = Min[xs]; xmax = Max[xs]; ymin = Min[ys]; ymax = Max[ys];

gbase = RegionPlot[L[1, x, y] > 0 && L[2, x, y] > 0 && L[3, x, y] > 0, {x, xmin, xmax},{y, ymin, ymax}]

gb=Table[Plot3D[L[i, x, y], {x, xmin, xmax}, {y, ymin, ymax}, RegionFunction -> Function[{x, y, z},

L[1, x, y] > 0 && L[2, x, y] > 0 && L[3, x, y] > 0]], {i, 3}]

gc=Table[ContourPlot[L[i, x, y], {x, xmin, xmax}, {y, ymin, ymax}, RegionFunction -> Function[{x, y, z},

L[1, x, y] > 0 && L[2, x, y] > 0 && L[3, x, y] > 0]], {i, 3}]

図

3:

三角形要素

かつらだ桂田

まさし

(17)

Mathematica のコード例と実行結果

授業

WWW

サイトからブラウザで、あるいは次のようにして

20230523.nb

を入手して実行してみよう。

ターミナルで入手して開く

curl -O https://m-katsurada.sakura.ne.jp/ana2023/20230523.nb open 20230523.nb

(18)

Mathematica のコード例と実行結果

図

4:

左から

L

0

, L

1

, L

2の鳥瞰図と等高線

かつらだ桂田

まさし

(19)

5.2.3 面積座標の積の積分と (L j , L i ) e

面積座標の積の積分については、便利な公式がある。

0

以上の任意の整数

i, j , k

に対して

(9)

Z Z

e

L 0 (x, y ) i L 1 (x, y) j L 2 (x, y) k dx dy = 2 | e | i!j!k!

(i + j + k + 2)! .

(20)

5.2.3 面積座標の積の積分と (L j , L i ) e

証明

P

0

= (0, 0), P

1

= (1, 0), P

2

= (0, 1)

で囲まれる三角形を

∆

とし、1次関数

φ: R

²

→ R

² を

φ(P

0

) = N

0

, φ(P

1

) = N

1

, φ(P

2

) = N

2

で定める。このとき

φ(∆) = e, φ(u, v) =

x

1

− x

0

x

2

− x

0

y

₁

− y

₀

y

₂

− y

₀

u v

+ x

0

y

₀

,

det φ

^′

(u, v) = det

x

₁

− x

₀

y

₁

− y

₀

x

₂

− x

₀

y

₂

− y

₀

= 2 | e | .

ゆえに変数変換

(x, y) = φ(u, v)

を行なうと

Z Z

e

L

ⁱ₀

L

^j₁

L

^k₂

dx dy = Z Z

∆

L

₀

(φ(u, v ))

ⁱ

L

₁

(φ(u, v ))

^j

L

₂

(φ(u, v ))

^k

| det φ

^′

(u, v) | du dv .

= 2 | e | Z Z

∆

L

₀

(φ(u, v ))

ⁱ

L

₁

(φ(u, v ))

^j

L

₂

(φ(u, v))

^k

du dv .

かつらだ桂田

まさし

(21)

5.2.3 面積座標の積の積分と (L j , L i ) e

L

i

(φ(u, v ))

^は

(1

^次関数と

1

次関数の合成であるから

)

^、

u, v

^{についての}

1

^{次関数で、}

L

i

(N

j

) = L

i

(φ(P

j

)) = δ

ij

を満たすことから

L

0

(φ(u, v )) = 1 − u − v , L

1

(φ(u, v )) = u, L

2

(φ(u, v )) = v

(

各等式の両辺は

1

次関数で、

N

j での値が一致するから、全体で一致する

).

ゆえに

Z Z

e

L

ⁱ₀

L

^j₁

L

^k₂

dxdy = 2 | e | Z Z

∆

(1 − u − v )

ⁱ

u

^j

v

^k

dudv

= 2 | e | Z

1

0

u

^j

Z

_1−u

0

(1 − u − v )

ⁱ

v

^k

dv

du.

右辺の内側の積分で、

v = (1 − u)t (0 ≤ t ≤ 1)

と変数変換すると

dv = (1 − u)dt, (1 − u − v )

ⁱ

= ((1 − u) − (1 − u)t)

ⁱ

= (1 − u)

ⁱ

(1 − t)

ⁱ であるから

Z

_1−u

0

(1 − u − v)

ⁱ

v

^k

dv = Z

1

0

(1 − u)

ⁱ

(1 − t)

ⁱ

(1 − u)

^k

t

^k

· (1 − u)dt

= (1 − u)

^i+k+1

Z

₁

0

(1 − t)

ⁱ

t

^k

dt.

(22)

5.2.3 面積座標の積の積分と (L j , L i ) e

ゆえに

Z Z

e

L

ⁱ₀

L

^j₁

L

^k₂

dx dy = 2 | e | Z

₁

0

(1 − u)

^i+k+1

u

^j

du Z

₁

0

(1 − t)

ⁱ

t

^k

dt

= 2|e|B(i + k + 2, j + 1)B (i + 1, k + 1).

ただし

B

は次式で定義されるベータ関数である

: B(p, q) =

Z

1 0

(1 − t)

^p−1

t

^q−1

dt (p, q > 0).

このベータ関数と、ガンマ関数

Γ (x ) :=

Z

_∞

0

t

^x⁻¹

e

⁻^t

dt (x > 0)

について、次の有名な公式が成り立つ

(

証明は例えば桂田

[?]

§

E.2

の命題

E.2.4)

。

B(p, q) = Γ (p)Γ (q)

Γ (p + q) , Γ (n + 1) = n! (n ∈ N ).

ゆえに

Z Z

e

L

ⁱ0

L

^j₁

L

^k2

dx dy = 2|e| Γ (i + k + 2) Γ (j + 1)

Γ (i + j + k + 3) · Γ (i + 1)Γ (k + 1) Γ (i + k + 2)

= 2|e| Γ (i + 1)Γ (j + 1)Γ (k + 1)

Γ (i + j + k + 3) = 2|e| i !j!k!

(i + j + k + 2)! .

かつらだ桂田

まさし

(23)

5.2.3 面積座標の積の積分と (L j , L i ) e

i = j

の場合、それ以外の添字

( ∈ { 0, 1, 2 } )

を

k, ℓ

とすると

(L

j

, L

i

)

e

= Z Z

e

L

²_i

L

⁰_k

L

⁰_ℓ

dxdy = 2 | e | 2!0!0!

(2 + 0 + 0 + 2)! = 1 6 | e | . i ̸ = j

の場合、それ以外の添字

( ∈ { 0, 1, 2 } )

を

k

とすると

(L

_j

, L

_i

)

_e

= Z Z

e

L

¹_i

L

¹_j

L

⁰_k

dxdy = 2 | e | 1!1!0!

(1 + 1 + 0 + 2)! = 1 12 | e | .

(24)

5.2.4 L i (x , y ) = a i + b i x + c i y の係数決定と ⟨ L j , L i ⟩ e

L

i

(x , y ) = a

i

+ b

i

x + c

i

y

とおくと

δ

ij

= L

j

(x

i

, y

i

) = a

j

+ b

j

x

i

+ c

j

y

i

= (1 x

i

y

i

)



 a

j

b

j

c

j





であるから

(♡)



 1 0 0 0 1 0 0 0 1



 =



 1 x

0

y

0

1 x

1

y

1

1 x

2

y

2







 a

0

a

1

a

2

b

0

b

1

b

2

c

0

c

1

c

2



 .

A :=



 1 x

0

y

0

1 x

1

y

1

1 x

2

y

2





とおくと

det A = det



 1 x

0

y

0

0 x

1

− x

0

y

1

− y

0

0 x

2

− x

0

y

2

− y

0



 =

x

1

− x

0

y

1

− y

0

x

2

− x

0

y

2

− y

0

= 2 | e | > 0.

ゆえに

A

は逆行列を持つ。

かつらだ桂田

まさし

(25)

5.2.4 L i (x , y ) = a i + b i x + c i y の係数決定と ⟨ L j , L i ⟩ e

(♡)

^から



 a

0

a

1

a

2

b

0

b

1

b

2

c

0

c

1

c

2



 = A

⁻¹

.

Cramer

の公式によって

A

⁻¹

= 1 det A



 

 

( − 1)

¹⁺¹

x

1

y

1

x

2

y

2

( − 1)

¹⁺²

1 y

1

1 y

2

( − 1)

¹⁺³

1 x

1

1 x

2

(−1)

²⁺¹

x

0

y

0

x

2

y

2

(−1)

²⁺²

1 y

0

1 y

2

(−1)

³⁺³

1 x

0

1 x

2

( − 1)

³⁺¹

x

0

y

0

x

1

y

1

( − 1)

³⁺²

1 y

0

1 y

1

( − 1)

³⁺³

1 x

0

1 x

1



 

 

⊤

= 1

det A



 x

1

y

2

− y

1

x

2

−(y

2

− y

1

) x

2

− x

1

− (x

0

y

2

− y

0

x

2

) y

2

− y

0

− (x

2

− x

0

) x

0

y

1

− y

0

x

1

− (y

1

− y

0

) x

1

− x

0





⊤

= 1

det A



 x

1

y

2

− y

1

x

2

x

2

y

0

− y

2

x

0

x

0

y

1

− y

0

x

1

y

1

− y

2

y

2

− y

0

y

0

− y

1

x

2

− x

1

x

0

− x

2

x

1

− x

0



 .

(26)

5.2.4 L i (x , y ) = a i + b i x + c i y の係数決定と ⟨ L j , L i ⟩ e

A

⁻¹の下

2

行

b

k

, c

k

(k = 0, 1, 2)

が得られれば、

⟨ L

j

, L

i

⟩

_e が計算できる:

⟨ L

j

, L

i

⟩

_e

= Z Z

e

∇ L

j

· ∇ L

i

dx dy = Z Z

e

b

j

c

j

· b

i

c

i

dx dy (10)

= (b

j

b

i

+ c

j

c

i

) | e | .

かつらだ桂田

まさし

(27)

5.3 要素係数行列の計算

「積分は積分範囲を分割して計算し、後から和を取ればよい」ので、弱形式

⟨b u, v b ⟩ = (f , v b ) + [g

2

, v b ] (b v ∈ X b )

は

(11)

N_e

X

k=1

⟨b u, v b ⟩

_e_k

=

N_e

X

k=1

(f , b v )

e_k

+ [g

2

, b v] ( b v ∈ X b )

となる。ただし

⟨b u, b v ⟩

_e_k

:=

Z

e_k

∇b u(x ) · ∇b v (x ) dx , (f , v b )

e_k

:=

Z

e_k

f (x) b v(x) dx, [g

2

, v b ] :=

Z

Γ₂

g

2

v dσ. b

そこで

u

^j

:= u(N b

j

), v

^j

:= b v (N

j

) (j = 0, 1, 2)

を用いて、

b u =

X

2 j=0

u

^j

L

j

, v b = X

2

i=0

v

ⁱ

L

i

(e

k 上

)

と表すと

(12) ⟨b u, b v ⟩

_e_k

= X

2

j=0

X

2 i=0

v

ⁱ

A

^(k)_ij

u

^j

, (f , v b )

e_k

= X

2

i=0

v

ⁱ

f

_i^(k)

,

ただし

(13) A

^(k)_ij

:= ⟨ L

j

, L

i

⟩

e_k

, f

_i^(k)

:= (f , L

i

)

e_k

.

(28)

5.3 要素係数行列の計算

そこで

(14a) u

k

:=



 u

⁰

u

¹

u

²



 , v

k

:=



 v

⁰

v

¹

v

²



 , f

k

:=



  f

₀^(k⁾

f

₁^(k⁾

f

₂^(k⁾



  ,

(14b) A

_k

:=



 

A

^(k)₀₀

A

^(k)₀₁

A

^(k)₀₂

A

^(k)₁₀

A

^(k)₁₁

A

^(k)₁₂

A

^(k)₂₀

A

^(k)₂₁

A

^(k)₂₂



 

とおくと、

A

k は対称行列で、

(15) ⟨b u, v b ⟩

e_k

= v

_k^T

A

k

u

k

, (f , b v)

e_k

= v

_k^T

f

k

.

(

線積分

[g

2

, v] ˆ

については、今回の授業では説明を省略する。とりあえず

g

2

= 0

と考えて授業を聴いて下さい。桂田

[?]

には書いてある。

)

A

^(k)_ij

= ⟨ L

_j

, L

_i

⟩

_e_k の計算は

(10)

で済んでいる。

かつらだ桂田

まさし

(29)

5.3 要素係数行列の計算具体的な成分の計算

f i (k) = (f , L i ) e

_k については、例えば

L i

を

1

次関数補間して

(f , L i ) e

_k

≒



 X 2

j=0

f (N j )L j , L i





e

k

= X 2

j=0

f (N j ) (L j , L i ) e

k

.

のように近似すれば

(L j , L i ) e

k

=

| e k | /6 (i = j

のとき

),

| e k | /12 (i ̸ = j

^のとき

).

であるから

f 0 (k) ≒ | e k |

12 (2f (N 0 ) + f (N 1 ) + f (N 2 )), f 1 (k) ≒ | e k |

12 (f (N 0 ) + 2f (N 1 ) + f (N 2 )), f 2 (k) ≒ | e k |

12 (f (N 0 ) + f (N 1 ) + 2f (N 2 )).

(30)

5.4 近似方程式の組み立て — 直接剛性法

u

i

:= u(P b

i

), v

i

:= b v (P

i

) (i = 1, 2, · · · , m)

^として、

(16) u :=



 

  u

1

u

2

. . . u

m



 

  , v :=



 

  v

1

v

2

. . . v

m



 

 

とおく。

有限要素

e

k

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

e

)

の節点

N

0

, N

1

, N

2 に対して、

N

0

= P

i_k,0

, N

1

= P

i_k,1

, N

2

= P

i_k,2

となる整数

i

k,0

, i

k,1

, i

k,2を取る

(

これらを全体節点番号と呼ぶ

)

。

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)

^⊤

.

このように

f

_k^∗ を定義すると、次が成り立つ。

(17) (f , v b )

e_k

= v

_k^⊤

f

k

= v

^⊤

f

_k^∗

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

e

).

かつらだ桂田

まさし

(31)

5.4 近似方程式の組み立て — 直接剛性法

同様の考え方で、行列

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

で定める。例えば

i

k,0

< i

k,1

< i

k,2ならば

A

^∗k

=



 

 

i

k,0

i

k,1

i

k,2

↓ ↓ ↓

A

^(k)₀₀

A

^(k)₀₁

A

^(k)₀₂

← i

k,0

A

^(k)₁₀

A

^(k)₁₁

A

^(k)₁₂

← i

k,1

A

^(k)₂₀

A

^(k)₂₁

A

^(k)₂₂

← i

k,2



 

 

(

書いてない成分は

0).

これを用いると

(18) ⟨ u, b v b ⟩

e_k

= v

k^⊤

A

k

u

k

= v

^⊤

A

^∗k

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

e

).

(32)

5.4 近似方程式の組み立て — 直接剛性法

ゆえに弱形式

(11)

は

(g

2

= 0

と考えている

)

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

^∗

:=

N_e

X

k=1

A

^∗_k

, f

^∗

:=

N_e

X

k=1

f

_k^∗

とおけば

(19) v

^⊤

(A

^∗

u − f

^∗

) = 0.

かつらだ桂田

まさし

(33)

5.4 近似方程式の組み立て — 直接剛性法

ここで

v

は

Y := {



  v

1

. . . v

m



  ∈ R

^m

; v

i

= 0 (P

i

∈ Γ

1なる

i ) }

の任意の元であるから、

(19)

^{は次と同値である。}

A

^∗

u − f ∈ Y

^⊥

= {



  w

1

. . . w

m



  ∈ R

^m

; w

i

= 0 (P

i

̸∈ Γ

1なる

i)}

すなわち

(20) A

^∗∗

u = f

^∗∗

.

ここで

A

^∗∗

:= A

^∗の第

i

行

(P

i

∈ Γ

1なる

i )

を除いた行列

, f

^∗∗

:= f

^∗の第

i

成分

(P

i

∈ Γ

1なる

i )

を除いた縦ベクトル

.

(34)

5.4 近似方程式の組み立て — 直接剛性法

I := { i ∈ N | 1 ≤ i ≤ m } , I

1

:= { i ∈ I | P

i

∈ Γ

1

}

とおく

(I

₁ は

Γ

₁上にある節点の節点番号全体の集合)。

条件

u

_i

= g

₁

(P

_i

) (i ∈ I

₁

)

があるから、これを代入して

u

i

(i ∈ I

1

)

を消去できる。以下それを実行する。

A

^∗∗ を列ベクトルで

A

^∗∗

= (a

1

a

2

· · · a

m

)

のように表示すると、(20)は

(21)

X

m i=1

a

i

u

_i

= f

^∗∗

.

左辺の

X

m

i=1

を

X

i∈I\I₁

+ X

i∈I1

と分解して、移項すると

X

i∈I\I1

a

i

u

i

= f

^∗∗

− X

i∈I₁

a

i

u

i

.

かつらだ桂田

まさし

(35)

5.4 近似方程式の組み立て — 直接剛性法

ゆえに

Au ∗ = f ,

ただし

u ∗ := u

の第

i

成分

(i ∈ I 1 )

を除いた縦ベクトル

, A := A ∗∗

の第

i

列

(i ∈ I 1 )

を除いた正方行列

, f := f ∗∗ − X

i ∈ I

1

a i u i .

実際にプログラムを作成するとき、

A

^や

f

が容易に求められることは次回解説する。

(36)

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

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

1

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

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

Γ 1 = { (x, y) | x = 0, 0 ≤ y ≤ 1 } ∪ { (x, y) | 0 ≤ x ≤ 1, y = 0 } , Γ 2 = { (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 Γ 1 , ∂u

∂n = 0 on Γ 2 .

かつらだ桂田

まさし

(37)

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

- 6

O x

y

x = 1 y = 1

Ω

図

5:

領域

Ω

- 6

x y

P 0

P 1

P 2

P 3

P 4 P 5

P 6 P 7

P 8

e 0 e 1

e 2 e 3

e 4 e 5

e 6 e 7

図

6:

要素分割

(38)

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

正方形領域

Ω

を図

2

のように

3

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

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

(タイプ I, II

と呼ぶことにする)。各々に図

3

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

(左下から反時計回り)。

タイプ

I e

0

, e

2

, e

4

, e

6

.

タイプ

II e

1

, e

3

, e

5

, e

7

.

N

0

N

1

N

₂

Type I

N

0

N

₁

N

₂

Type II

図

7:

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

かつらだ桂田

まさし

(39)

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

タイプが同じならば、要素係数行列

A

_k

,

要素自由項ベクトル

f

k が等しいことはすぐ分かる。それぞれ

A

_I

, A

_II

, f

_I

, f

_IIで表すことにする。

タイプ

I

については

D = h

²

, S = h

²

2 ,

A

_I

= 1 2



 1 − 1 0

− 1 2 − 1 0 − 1 1



 , f

I

= 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



 , f

II

= f h

²

6



 1 1 1



 .

(40)

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

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

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

(

スライド見る場合も手で写すことを推奨

—

要素タイプは不要

)

。

要素

e

0

e

1

e

2

e

3

e

4

e

5

e

6

e

7

要素タイプ

I II I II I II I II N

0の全体節点番号

0 0 1 1 3 3 4 4 N

3 4 4 5 6 7 7 8 N

4 1 5 2 7 4 8 5

N

0

N

1

N

2

Type I

N

0

N

1

N

2

Type II

- 6

x y

P

0

P

1

P

2

P

3

P

4

P

5

P

6

P

7

P

8

e

0

e

1

e

2

e

3

e

4

e

5

e

6

e

7

かつらだ桂田

まさし

(41)

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

これから

Galerkin

^{法の弱形式は}

v

^⊤

1 2



 



2 − 1 0 − 1 0 0 0 0 0

− 1 4 − 1 0 − 2 0 0 0 0

0 −1 2 0 0 −1 0 0 0

− 1 0 0 4 − 2 0 − 1 0 0

0 −2 0 −2 8 −2 0 −2 0

0 0 − 1 0 − 2 4 0 0 − 1

0 0 0 −1 0 0 2 −1 0

0 0 0 0 − 2 0 − 1 4 − 1

0 0 0 0 0 −1 0 −1 2



 





 

 u

0

u

1

u

2

u

3

u

4

u

5

u

6

u

7

u

8



 



= v

^⊤

f h

²

6



 

 2 3 1 3 6 3 1 3 2



 



for

∀ v ∈ n

(v

0

, v

1

, v

2

, v

3

, v

4

, v

5

, v

6

, v

7

, v

8

)

^⊤

∈ R

⁹

v

0

= v

1

= v

2

= v

3

= v

6

= 0 o

.

(42)

5.5 連立 1 次方程式の具体例 Dirichlet 境界条件の考慮

P

i

∈ Γ

1となる

i

^について

(

^今の場合

i = 0, 1, 2, 3, 6)

^、第

i

^{行は削除してよい。}

1 2



 



0 − 2 0 − 2 8 − 2 0 − 2 0

0 0 − 1 0 − 2 4 0 0 − 1

0 0 0 0 −2 0 −1 4 −1

0 0 0 0 0 − 1 0 − 1 2



 





 

 u

0

u

1

u

2

u

3

u

4

u

5

u

6

u

7

u

8



 



= f h

²

6



 

 6 3 3 2



 

 .

また

P

i

∈ Γ

1となる

i

について、

u

i

= g

1

(P

i

) (= 0).

これを代入して移項すると

1

2



 



8 −2 −2 0

− 2 4 0 − 1

− 2 0 4 − 1

0 −1 −1 2



 





 

 u

4

u

5

u

7

u

8



 

 = f h

²

6



 

 6 3 3 2



 

 +



 



g

1

(P

1

) + g

1

(P

3

) g

1

(P

2

)/2 g

1

(P

6

)/2

0



 

 .

すなわち



 



4 − 1 − 1 0

−1 2 0 −1/2

− 1 0 2 − 1/2 0 − 1/2 − 1/2 1



 





 

 u

4

u

5

u

7

u

8



 

 =



 

 f h

²

f h

²

/2 f h

²

/2 f h

²

/3



 

 .

かつらだ桂田

まさし

(43)

5.5 連立 1 次方程式の具体例 Dirichlet 境界条件の考慮

上のやり方では、係数行列とベクトルが

“

縮小

”

される。いくつか留意すべき点

:

例えば

MATLAB

^では、

(0, 1, 2, 3, 6), (4, 5, 7, 8)

^{という添字ベクトル} を用意すれば、全体係数行列と全体自由項ベクトルを求めるのは

(

コーディング上は

)

容易である。

自分で疎行列を扱うコードを書いていたりする場合はそれなりに面倒。

データの移動にも計算コストがかかる。

Dirichlet

境界条件の処理には、他のやり方

(

行列とベクトルを縮小しな

い方法

)

^もある。

(44)

5.5 連立 1 次方程式の具体例 Dirichlet 境界条件の考慮

ベクトル、行列の縮小を避ける方法

(1)

P

i

∈ Γ

1となる

i (

この例では

i = 0, 1, 2, 3, 6)

に対して

i

番目の方程式

(v

^⊤のせいで

0

^⊤

u = 0)

を

u

i

= g

1

(P

i

)

で置き換える。

(

結果的に係数行列の第

i

行は

e

_i^⊤で置き換える

)

とすることで

v

^⊤が外せる。



 



1 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0

0 − 1 0 − 1 4 − 1 0 − 1 0

0 0 −1/2 0 −1 2 0 0 −1/2

0 0 0 0 0 0 1 0 0

0 0 0 0 −1 0 −1/2 2 −1/2

0 0 0 0 0 − 1/2 0 − 1/2 1



 





 

 u

0

u

1

u

2

u

3

u

4

u

5

u

6

u

7

u

8



 



=



 

 g

1

(P

0

) g

1

(P

1

) g

1

(P

2

) g

1

(P

3

) f h

²

f h

²

/2 g

1

(P

6

) f h

²

/2 f h

²

/3



 

 .

これは正しい方程式であるが、係数行列が対称でなくなっている

(

数値計算で不利

)

。そこで、係数行列の

i

列に

0

でない非対角要素があれば、それと

g

1

(P

i

)

との積を右辺

に移項する。

(

その結果は次のスライド

)

かつらだ桂田

まさし

(45)

5.5 連立 1 次方程式の具体例 Dirichlet 境界条件の考慮



 



1 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0

0 0 0 0 4 − 1 0 − 1 0

0 0 0 0 − 1 2 0 0 − 1/2

0 0 0 0 0 0 1 0 0

応用数値解析特論第 6