Astroの光線のサムネイル。

pubDate: 2024-03-18

author: sakakibara

有限要素法

行列表現

ベクトル解析

勾配・発散・回転

以下の記号を知っているだろうか。

gradϕdivvrotv\mathrm{grad} \phi \\ \mathrm{div} \mathbf{v} \\ \mathrm{rot} \mathbf{v}

これらは所謂ベクトル解析と呼ばれる分野に登場する記号だ。
grad\mathrm{grad}から説明しよう。

勾配

ϕ\phiをnn次元実数空間x∈Rn\mathbf{x}\in \mathbb{R}^nから実数R\mathbb{R}への関数とする。 形式的に記述すると以下のような感じだ。

ϕ:Rn→R\phi : \mathbb{R}^n \rightarrow \mathbb{R}

例としては等高線が書かれた地図のようなものや空間に分布する温度などが取り上げられる。 (この手の関数、 つまり、ベクトルから実数への関数はかなり多くの研究対象となっている。)

空間の温度の分布について考える際、その位置に対する温度の変位は指定した位置による微分で定義されるだろう。 位置を33次元実数空間x=(x1,x2,x3)\mathbf{x}=(x_1, x_2, x_3)で表すと、

(∂ϕ∂x1,∂ϕ∂x2,∂ϕ∂x3)\bigg(\frac{\partial{\phi}}{\partial{x_1}}, \frac{\partial{\phi}}{\partial{x_2}}, \frac{\partial{\phi}}{\partial{x_3}}\bigg)

となる。 または∇=(∂∂x1,∂∂x2,∂∂x3)\nabla = (\frac{\partial}{\partial{x_1}}, \frac{\partial}{\partial{x_2}}, \frac{\partial}{\partial x_3})を用いて

∇ϕ\nabla\phi

のように表される。

そこで、以下のように演算子を定義する。

grad:=∇\mathrm{grad} := \nabla

発散

v\mathbf{v}をnn次元実数空間x∈Rn\mathbf{x}\in \mathbb{R}^nからベクトル空間VVへの関数とする。 形式的に記述すると以下のような感じだ。

v:Rn→V\mathbf{v} : \mathbb{R}^n \rightarrow V

例としては水流の分布などを想像してほしい。

水流の湧き出しについて考える際、湧き出しの量はその位置の各方向で微小面積あたりどれほど水量が変化(水がその方向から出ていって入った量)したかで定義されるだろう。 位置を33次元実数空間x=(x1,x2,x3)\mathbf{x}=(x_1, x_2, x_3)で表すと、

∂v1∂x1+∂v2∂x2+∂v3∂x3\frac{\partial{v_1}}{\partial{x_1}} + \frac{\partial{v_2}}{\partial{x_2}} + \frac{\partial{v_3}}{\partial{x_3}}

∇\nablaを用いて表すと

∇⋅v\nabla \cdot \mathbf{v}

のように表される。

そこで、以下のように演算子を定義する。

div:=(∇⋅)\mathrm{div} := (\nabla \cdot)

回転

水流の回転について考える際、回転の量はその位置の各方向の変位に対してどれほど水量がその方向と垂直方向に変化したかで定義されるだろう。 位置を33次元実数空間x=(x1,x2,x3)\mathbf{x}=(x_1, x_2, x_3)で表すと、

(∂v3∂x2−∂v2∂x3,∂v1∂x3−∂v3∂x1,∂v2∂x1−∂v1∂x2)\bigg(\frac{\partial{v_3}}{\partial{x_2}} - \frac{\partial{v_2}}{\partial{x_3}}, \frac{\partial{v_1}}{\partial{x_3}} - \frac{\partial{v_3}}{\partial{x_1}}, \frac{\partial{v_2}}{\partial{x_1}} - \frac{\partial{v_1}}{\partial{x_2}}\bigg)

∇\nablaを用いて表すと

∇×v\nabla \times \mathbf{v}

のように表される。

そこで、以下のように演算子を定義する。

rot:=(∇×)\mathrm{rot} := (\nabla \times)

それぞれの関係式

div\mathrm{div}, rot\mathrm{rot}, grad\mathrm{grad}それぞれを結ぶ重要な関係式がある。 それが以下だ。

rot(grad)=0div(rot)=0\begin{aligned} \mathrm{rot}(\mathrm{grad}) &= 0 \\ \mathrm{div}(\mathrm{rot}) &= 0 \end{aligned}

すなわち回転するような勾配は存在しないし、発散するような回転は存在しないということである。
これは実際に∇\nablaを計算してみるとわかる。 また、以下の演算をラプラシアンやラプラス演算子などと呼ぶ

Δ=div(grad)\Delta = \mathrm{div}(\mathrm{grad})

これは調和作用素などと呼ばれることもあり、物理の文脈でよく出てくる。
また、なぜだかこの逆には特別な名前がついていない。

?=grad(div)? = \mathrm{grad}(\mathrm{div})

この演算には特別な名前ついていないものの、以下の様に重要な役割を果たす。

rot(rot)=[grad,div]=grad(div)−div(grad)=grad(div)−Δ\begin{aligned} \mathrm{rot}(\mathrm{rot}) &= [\mathrm{grad}, \mathrm{div}] \\ &= \mathrm{grad}(\mathrm{div}) - \mathrm{div}(\mathrm{grad}) \\ &= \mathrm{grad}(\mathrm{div}) - \Delta \end{aligned}

この演算には名前が無いので個人的にロトロトと読んでいる。

ベクトル場と離散表現

さて、これらの演算子はベクトル場やスカラー場に対して定義されるものである。 以下ではこれらの演算子を行列を用いて表現する方法について考える。 以下の内容はA.Bossavit, “Computational Electromagnetism: Variational Formulations, Complementarity, Edge Elements”に基づいている。

勾配行列

さて、以上でベクトル場v\mathbf{v}やスカラー場ϕ\phiとそれにまつわる演算子について軽く触れた。

先程触れた勾配・発散・回転と同様の概念はネットワーク的に表現できる。 ネットワーク的にとは接続行列や隣接行列を使用して定義できるのだ。

頂点n1n_1から辺eeが出て、頂点n2n_2へ入っているとする。 この関係を

∂e=n2−n1\partial{e} = n_2 - n_1

のように表現する。 これを様々な辺について考える。
辺eke_kが頂点nin_iから出て頂点njn_jへ出ているとする。

∂ek=∑iGkini\partial{e_k} = \sum_i G_{ki} n_i

ここでGkiG_{ki}は

Gki={辺ekが頂点niへ入っているのであれば1辺ekが頂点niから出ているのであれば−1other0G_{ki} = \begin{cases} \text{辺}e_k \text{が頂点} n_i \text{へ入っているのであれば}1 \\ \text{辺}e_k \text{が頂点} n_i \text{から出ているのであれば}-1 \\ \mathrm{other} 0 \end{cases}

のようにして定義される。

このような要素を持つ行列をG=[Gki]\mathrm{G}=[G_{ki}]と定義する。
これを勾配行列という。

頂点の個数をNN, 辺の個数をEEとすると勾配行列G\mathrm{G}はE×NE \times Nの行列である。 頂点に値を振り宛てたNN個の要素をもつベクトルをϕ∈RN\phi \in \mathbb{R}^Nとすると、 その勾配はGϕ\mathrm{G}\phiのようにして表される。

勾配行列G\mathrm{G}は所謂接続行列になっている。 そして勾配行列の積をdivgrad\mathrm{div}\mathrm{grad}とのアナロジーからラプラス行列という。

L=GTG\mathrm{L} = \mathrm{G}^T\mathrm{G}

回転行列

勾配行列で考えた辺と頂点に対する考えをさらに拡張し、 面と辺に対しても同じような考えを適用してみる。

面ffの境界の辺e1,e2,e3,e4e_1, e_2, e_3, e_4が面を囲む向きに定まっているとする。この関係を

∂f=e1+e2+e3+e4\partial{f} = e_1 + e_2 + e_3 + e_4

のように表現する。

これを様々な面について考える。
面fkf_kが辺eie_iに囲まれているとする。

∂fk=∑iRkiei\partial{f_k} = \sum_i R_{ki} e_i

ここでRkiR_{ki}は

Rki={面fkが辺eiによって適切な向きによって囲まれているのであれば1面fkが辺eiによって適切な向きの逆向きによって囲まれているのであれば−1other0R_{ki} = \begin{cases} \text{面}f_k\text{が辺}e_i\text{によって適切な向きによって囲まれているのであれば}1 \\ \text{面}f_k\text{が辺}e_i\text{によって適切な向きの逆向きによって囲まれているのであれば}-1 \\ \mathrm{other} 0 \end{cases}

のようにして定義される。

このような要素を持つ行列をR=[Rki]\mathrm{R}=[R_{ki}]と定義する。
これを回転行列という。

辺の個数をEE, 面の個数をFFとすると回転行列R\mathrm{R}はF×EF \times Eの行列である。 辺に値を振り宛てたEE個の要素をもつベクトルをv∈RE\mathbf{v} \in \mathbb{R}^Eとすると、 その回転はRv\mathrm{R}\mathbf{v}のようにして表される。

ラプラス行列のようにロトロト行列RTR\mathrm{R}^T\mathrm{R}について考えることもできる。

発散行列

発散に関しても以上と同じ関係が導かれる。 それをD\mathbf{D}とする。

勾配・回転・発散 に関する関係式

さて、興味深いことにベクトル解析で見られた勾配・回転・発散の関係が離散化した状態でも見られる。

つまり、

rot(grad)=0div(rot)=0\mathrm{rot}(\mathrm{grad}) = 0 \\ \mathrm{div}(\mathrm{rot}) = 0

が表現できる。

面fについて

∂f=e1+e2+e3+e4∂e1=n2−n1∂e2=n3−n2∂e3=n4−n3∂e4=n1−n4\begin{aligned} \partial{f} &= e_1 + e_2 + e_3 + e_4 \\ \partial{e_1} &= n_2 - n_1 \\ \partial{e_2} &= n_3 - n_2 \\ \partial{e_3} &= n_4 - n_3 \\ \partial{e_4} &= n_1 - n_4 \end{aligned}

のように表現できるとしよう。

さて、ここで以下のような考えをしてみる。

∂∂f=∂(e1+e2+e3+e4)=∂e1+∂e2+∂e3+∂e4=n2−n1+n3−n2+n4−n3+n1−n4=0\begin{aligned} \partial{\partial{f}} &= \partial{(e_1 + e_2 + e_3 + e_4)} \\ &= \partial{e_1} + \partial{e_2} + \partial{e_3} + \partial{e_4} \\ &= n_2 - n_1 \\ &+ n_3 - n_2 \\ &+ n_4 - n_3 \\ &+ n_1 - n_4 \\ &= 0 \end{aligned}

なんと、∂f\partial{f}の境界を取る作用をすると、それは00になる。

そしてこれを行列的な表現をすると、

RG=O\mathrm{RG} = O

となる。 また、発散行列D\mathrm{D}と回転行列R\mathrm{R}についても同じことが言える。

DR=O\mathrm{DR} = O

辺要素のrot

有限要素法では解析対象を要素と呼ばれる小さな領域に分割する。 解析対象に関する物理場は要素内での補間により表現する。 要素内での補間の方法には大きく３つある。

ベクトル・ポテンシャルなどを扱う際には辺補間を使う。
ベクトル・ポテンシャルA\mathbf{A}を辺補間で補間する際に磁束密度B\mathbf{B}は

B=rotA\mathbf{B} = \mathrm{rot}\mathbf{A}

のように表される。 そしてその辺要素の値をベクトルとした{A}\{A\}と{B}\{B\}について

{B}=R{A}\{B\} = R \{A\}

として表現できる。

B={∇×N}T{A}\mathbf{B} = \{{\nabla \times \mathbf{N}}\}^T\{A\}

また、磁束密度B\mathbf{B}は以下のようも表現できるとする。

B={M}T{B}={M}TR{A}={RT{M}}T{A}\begin{aligned} \mathbf{B} &= \{M\}^T\{B\} \\ &= \{M\}^TR\{A\} \\ &= \{R^T\{M\}\}^T\{A\} \end{aligned}

すなわち、

{∇×N}=RT{M}\{\nabla \times \mathbf{N}\} = R^T\{M\}

が成り立つ。

注意しなければならないこととして、面から辺は計算できるが、辺から面は計算できないことだ。 つまり、辺の形状関数が判明していれば点の形状関数が求められるが逆は必ずしも成り立たないということである。