計算数理工学論文集 Vol. 13 (2013 年 11 月), 論文 No. 17-131129

JASCOME

Maxwell 方程式に対する PMCHWT 定式化におけるHdiv 内積を用いた離散化について
On discretisation methods with Hdiv scalar productfor PMCHWT formulations for Maxwell’s equations
新納 和樹 1) ，西村 直志 2)
Kazuki NIINO and Naoshi NISHIMURA
1) 京都大学情報学研究科(〒 606-8501

京都市左京区吉田本町,

E-mail: niino@acs.i.kyoto-u.ac.jp)

2) 京都大学情報学研究科(〒 606-8501

京都市左京区吉田本町,

E-mail: nchml@i.kyoto-u.ac.jp)

A new discretisation method with Hdiv scalar product for boundary element methods
is proposed. We show that this method improves the accuracy of boundary element
methods for PMCHWT formulations for Maxwell’s equations. In particular, this method
is eﬀective for low frequency problems in which conventional methods suﬀer from badaccuracy. We verify the eﬃciency of this method through some numerical examples.
Key Words : Maxwell Equations, PMCHWT Formulation, Galerkin Method
1. 序論

こで本稿では Hdiv 内積を用いた Galerkin 法によって積分方(1)(2)といった程式を離散化することで, 周波数の値が小さい問題でも精度周期光学材料がストップバンド性や自然界に存在しない負のよく解が求まる数値解法を提案する. この解法では従来法で屈折率などの興味深い性質を実現できるとして盛んに研究さ計算する項に加え, 追加の項を計算する必要があるが, このれている. これらの光学材料は波長に対して微小な構造を数項は従来法で計算している値を再利用することで容易に求ま多く配置するため, そのような電磁波動散乱問題の解を数値るため, 増加する計算量は非常に小さいことがわかる.

近年, フォトニック結晶やメタマテリアル

本稿の構成は以下のとおりである. まず 2 節で問題とそれ的に高速, 高精度に求めるために境界要素法は有力な解法でを解くための積分方程式の定式化を行う. 次に 3 節で一般的ある.
境界要素法を用いた電磁波動散乱問題の解析に関しては,に使われている積分方程式の離散化法について述べ, 4 節で

既 に 多く の 研 究が行 われ て お り, 様々な定式化が 提案さ れ新しい Hdiv 内積を用いる離散化法を提案し, その利点を示ている. 特に波動散乱問題に対しては Poggio-Miller-Chang-す. その後 5 節で数値計算例を示し, 6 節で結論と今後の課

Harrington-Wu-Tsai(3) (PMCHWT) 定式化による解法が, 見題を述べる.

かけの固有値の影響を受けずに安定して解が得られる数値解法として知られている. しかしながら電磁波動散乱問題を様々なパラメータで安定して高速に解ける解法は, 現状ではまだ開発されていないと言える.
電磁波動散乱問題を様々なパラメータで解くことが難しい

2. 定式化
Fig. 1 に示す様に, 滑らかな境界 Γ に囲まれた領域を Ω− ,Ω+= R3 \Ω− とする. 考える問題は, 領域 Ω± 内で Maxwell方程式

原因の一つが低周波問題である. これは周波数 ω の零極限で∇ × E = iωµ± H,Eelectric Field Integral Operator (EFIO) が特異性を持つこ∇ × H = −iωε± Eとに起因する問題であり, PMCHWT 定式化においてもこれが原因で周波数が小さいときに解の精度が著しく悪化することが知られている. 低周波問題の解決法はいくつか提案されているが(4, 5, 6), いずれの方法も実装が複雑であったり, 問題の自由度が増加するなど計算量の増加を招いてしまう. そを満足し, 境界 Γ 上で境界条件

m := E + × n = E − × n,
j := n × H + = n × H −を満たす E, H を求めることである. ここに E は電場, H は2013 年 9 月 17 日受付，2013 年 10 月 22 日受理

磁場, ω は周波数, ε± , µ± はそれぞれ領域 Ω± における誘電nΩ次の式で定義される関数空間を考える:+Hdiv (Γ) = {t | t ∈ L2T (Γ), divS t ∈ L2 (Γ)}
Hcurl (Γ) = {t | t ∈ L2T (Γ), curlS t ∈ L2 (Γ)}.

Ω−

ここに L2T (Γ) は Γ 上で定義される接ベクトル場 t で∫(t, t)L2 (Γ) :=|t|2 dS < ∞x2TΓ
x3Γを満たす関数の集合,x1divS ϕ := −(∇ × (ϕ × n)) · ncurlS ϕ := ∇ · (ϕ × n)

Fig. 1

Boundary value problems.
である. 基底関数 ti (i = 1, · · · , N ) を Hdiv の要素とすると,率と透磁率, n は境界 Γ 上で定義される外向き単位本線ベク±トル, E , H±±はそれぞれ電場と磁場の Ω から Γ 上への極

Hdiv , Hcurl の定義より明らかに n × ti ∈ Hcurl である. 式 (1),(2) の両辺は Hdiv (Γ) に属すこと, Hdiv (Γ) と Hcurl (Γ) は双対空間であること (7) に注意すると, 式 (1), (2) を以下の様に限値とする. さらに Ω+ においてtesting するのが自然である.
E = E inc + E sca ,(n × ti , (P + + P − )m − iω(µ+ Q+ + µ− Q− )j)L2 (Γ)

H = H inc + H sca

Tとして散乱波 E sca , H sca に遠方での放射条件を課す. ここに E inc , H inc は入射波である.
する.

Ψkl =

(ekjl ∂j G(x),

Φkl =

j=1

Qil (x, y) =)L2 (Γ)

(4)

TT

=(ti , E inc )L2 (Γ)(5)

T(ti , iω(ε+ Φ+ + ε− Φ− )m + (Ψ+ + Ψ− )j)L2 (Γ)
T= − (ti , Heijk nj (x)Ψkl (x − y),j=1 k=13 ∑3∑= − (n × ti , n × Hinc

(ti , (Ψ+ + Ψ− )m − iω(µ+ Φ+ + µ− Φ− )j)L2 (Γ)G(x) は Helmholtz 方程式の基本解である. さらに核関数
Pil (x, y) =(n × ti , iω(ε+ Q+ + ε− Q− )m + (P + + P − )j)L2 (Γ)

さらに式 (3), (4) は)1
δkl + 2 ∂k ∂l G(x).
k

ここに ekjl は三次元の交代記号, δkj は Kronecker のデルタ,3 ∑3∑(3)TTこの問題の積分方程式を表すために以下の積分核を導入

3∑= − (n × ti , E inc × n)L2 (Γ)inc)L2 (Γ)(6)Tと書き換えられる. 未知関数 m, j をeijk nj (x)Φkl (x − y),m=j=1 k=1N∑
i=1m i ti ,j=N∑

j i ti ,i=1を定義し, Ψ, Φ, P, Q の四つの核関数に対して, Γ 上での積分作用素を∫P (x, y) · t(y) dSy(P t)(x) =と展開し, 式 (5), (6) に代入することで線型方程式が得られる.

x ∈ Γ,

Γのように同じ記号によって定義する. PMCHWT 定式化による境界積分方程式はこれらの記号を用いて以下の様に表される.

4. 新しい離散化法
本節では本稿で新しく導入する離散化法について述べる.

4.1. 概要(P + + P − )m − iω(µ+ Q+ + µ− Q− )j = −E inc × n(1)iω(ε+ Q+ + ε− Q− )m + (P + + P − )j = −n × H inc(2)

3 節で述べた従来的な離散化法では, Hdiv と Hcurl が L2T (Γ)
内積によって双対空間となることを用い, Galerkin 法を行った. 本稿で提案する新しい離散化法では, 関数空間 Hdiv が4.2 節で定義する内積によって Hilbert 空間となることを利用し, その内積による Galerkin 法を定式化する. これによっ

3. 従来の Galerkin 法による境界積分方程式の離散化
本節では Galerkin 法による式 (1), (2) の離散化において広く用いられている方法について述べる.て得られる線形方程式は従来法で得られる方程式には含まれない境界積分方程式の法線成分が陽に現れる式となるため,特に低周波問題において精度が改善することが期待される.

以下, 4.2 節で Hdiv 内積を用いた新しい離散化法の具体的と書き換えることができる. この式において 3 節で述べた従な定式化について述べ, 4.3 節で実際に低周波問題に対して来法と同様に, 未知関数 m, j を基底 ti で展開することで線この手法が有効であることを示す.

形方程式が得られる.

4.2. 新しい離散化法の定式化

式 (9), (10) の両辺一項目はそれぞれ積分方程式 (1), (2) の関数空間 Hdiv (Γ) は内積接線成分に対応し, 定数 c が掛かっている式 (9), (10) の両辺二項目はそれぞれ式 (2), (1) の法線成分に対応している. し(u, v)Hdiv (Γ) := (u, v)L2 (Γ) + c(divS u, divS v)L2 (Γ)

たがって従来法では積分方程式 (1), (2) の接線成分のみ用いTによって Hilbert 空間となる. ここに c は正の定数である. 通ていたところを, 新しい離散化法では接線成分, 法線成分の常 (·, ·)Hdiv (Γ) は c = 1 によって定義されるが, ここでは数値両方を用いることになる. そのため新しい離散化法では従来計算上の理由から定数 c を導入する. 定数 c の値の決め方に法から新たに積分の計算を行う必要が無いので, 両離散化法

ついては 4.3 節で述べる. この内積を用いて式 (1), (2) をによる係数行列構成のための計算時間の差は僅かであると考えられる.

(si , (P + + P − )m − iω(µ+ Q+ + µ− Q− )j)Hdiv (Γ)

4.3. 新しい離散化法の利点= − (si , E inc × n)Hdiv (Γ)(7)(si , iω(ε+ Q+ + ε− Q− )m + (P + + P − )j)Hdiv (Γ)

本節では 4 節で導入した離散化法により, いわゆる “低周波問題” の精度が改善することを示す.

= − (si , n × H inc )Hdiv (Γ)(8)

式 (5), (6) から得られる線型方程式の係数行列において,最も絶対値の大きい項は対角項である. そこで各作用素からと testing することが考えられる. ここに si ∈ Hdiv (Γ) は(ti , n × sj )L2 (Γ)
Tを ij 成分とする行列が正則となるような基底であり, ti の双対基底と呼ばれる (8) . 基底 ti と si の選び方としては, 例え得られる係数行列の対角項のオーダーを計算すると( ∫)ti , ∇G × tj dS ∼ O(h3 )Γ()∫ti , ω Gtj dS ∼ O(ωh3 ) = O(ωh · h2 )(ば ti を Rao-Wilton-Glisson 基底 (9) (RWG 基底) としたとき

ti ,si を Buﬀa-Christiansen 基底 (10) (BC 基底) とすればよいが,1ωΓ∫)∇∇Gtj dS∼O

Γ

( )
( 2)
hh=Oω
ωh(11)(12)(13)となる. ここに h はメッシュの一辺の代表長である. 特に式後に見るようにこの限りではない.
式 (7), (8) において testing の基底に ti ではなく si を用い(11) における ∇ を接線微分と法線微分に分けると接線微分ているのは以下の理由による. 式 (7), (8) では式 (3), (4) との項が 0 となるため, この ∇ は実質法線微分であることに注比較し, testing に用いる関数から n× が脱落した形になって意すること. 式 (11), (12), (13) より, ωh ≪ 1 のとき式 (13)いる. したがって testing に用いる接ベクトル場が Γ 上の接の値が式 (11), (12) より大きくなる. 線形方程式 (5), (6) の平面において 90 度回転しているため, 式 (7), (8) の testing

左辺の係数行列は式 (11), (12), (13) の左辺の項より成っていの基底に ti ではなく, その双対基底 si を用いている. 式 (7),るため, ωh ≪ 1 のとき, 式 (13) の左辺の項に近い値をとる(8) はことがわかる. 一方で, 式 (13) 中の作用素∫∇∇Gt dS

divS P ± u = −k± n · Φ± u2

ΓdivS Q± u = −n · Ψ± u

は明らかに divS t = 0 を満たす任意の関数 t に作用させると0 となるため, この作用素を離散化して得られる式 (13) の左辺の行列は非正則行列である. したがって ωh ≪ 1 のとき, 式であることに注意すると(5), (6) を数値的に解くと精度が悪化することが予想される.
(si , (P + + P − )m − iω(µ+ Q+ + µ− Q− )j)L2 (Γ)

これが従来法に見られる低周波問題である.

T−−+c(divS si , n · {iω(µ Ψ + µ Ψ )j++

一方, 式 (9), (10) から得られる線型方程式の係数行列の対− (k+ Φ+ + k− Φ− )m})L2 (Γ)
22= − (si , E inc × n)L2 (Γ) − iωµ+ c(divS si , n · H inc )L2 (Γ)T(9)+−+−(si , iω(ε Q + ε Q )m + (P+−+ P )j)L2 (Γ)(T+c(divS si , n · {−(k + Φ+ + k− Φ− )j2−2si ,−− iω(ε Ψ + ε Ψ )m})L2 (Γ)++角項を計算すると, 両式の第一項目は)(∫∇G × tj dS ∼ O(h3 )si , n ×Γ)(∫Gtj dS ∼ O(ωh3 ) = O(ωh · h2 )si , ωn ×= − (si , n × H inc )L2 (Γ) − iωε+ c(divS si , n · E inc )L2 (Γ)T(10)

1n×ω∫Γ

∇∇Gtj dSΓ)
∼O( )
( 2)
hh=O
ωωh(14)(15)(16)

となり, 式 (11), (12), (13) と同じ値になるが, 第 2 項目は()∫c divS si , ωn ·∇G × tj dS ∼ O(cωh)
(17)
Γ()
∫
c divS si , ω 2 n ·Gtj dS ∼ O(cω 2 h2 )(18)
Γ()∫c divS si , n ·∇∇Gtj dS ∼ O(ch)(19)Γとなる. 特に式 (14) における ∇ と式 (19) における一つ目の∇ は式 (11) と同様, 実質法線微分であることに注意すること. 式 (14)∼(16) において式 (14) と式 (15) は ωh ≪ 1 のとき, 式 (16) と比較して微小項と見なすことができるため, 主要項は式 (16) であると言える. 同様に式 (17)∼(19) の中では式 (18) が式 (17) と ωh のオーダーで比較すると, 微小項と見求めた数値解の相対誤差(∫)∫
1|j c − j e |dS/H0 +|mc − me |dS/E02|Γ|ΓΓ
を様々な周波数に対してプロットした図である. ここに j c , mc
はそれぞれ j, m の数値解, j e , me は解析解, H0 = |H inc |,E0 = |E inc | である. 計算した周波数域の大部分において
c = c0 の Hdiv (Γ) 内積による解法が最も良い精度であることがわかる. 特に周波数が小さい問題において高い精度を示している. c = 1 の Hdiv (Γ) 内積による解法もある程度良い精度を示しているが, 周波数の値が小さい問題では, c = c0 の解法と比較して誤差が大きくなっている. 一方で従来法は特に周波数が小さい問題において精度が著しく悪化していないことがわかる.

なすことが出きるため, 式 (17) と式 (19) を主要項と見なす
Hdiv scalar product, c=c0
Hdiv scalar product, c=1

ことができる. これらについて

h2= cωh + ch
ωh1

(20)
0.1

題においても精度良く解ける線形方程式が構成できることが期待される. 式 (20) を c について解くと1
ω(ω + 1)

relative errorを満たすように定数 c を決めることで ωh ≪ 1 なる低周波問

c = c0 :=Conventional method

0.01
0.001
0.0001(21)

1e-05となる.
1e-06
0.00010.001

5. 数値計算例
一辺の長さ 1 の立方体形状の散乱体に平面波√iω ε+ µ+ v·x
E inc = E inc,0 e√iω ε+ µ+ v·x
H inc = H inc
0 e(22)

Fig. 20.01
frequency

0.1

1

Relative error for several frequencies withthe mesh of 384 quadrilateral elements.

(23)が入射する問題を考える. v = t (0, 0, 1), E inc= t (1, 0, 0),0

Hdiv scalar product, c=c0
Hdiv scalar product, c=1H inc= t (0, 1, 0) とする. ε± = µ± = 1 とすると解は入射波
01(22), (23) となるので, これを使って様々な値の周波数 ω に対して数値解の精度を確かめた.

0.1

れる c0 としたもの, 同じく c = 1 としたもの, 3 節で述べた従来的な Galerkin 法により得られる方程式 (5), (6) を解く数値解法について比較を行った. 全ての解法において散乱体表面

relative error

本稿では 4 節で導入した Hdiv (Γ) 内積によって得られる方程式 (9), (10) を解く解法において定数 c を式 (21) で定義さConventional method

0.01
0.001
0.0001
1e-05を四角形メッシュで分割し, その上で定義される roof top 関数を基底 ti として用いた. これによって si も roof top 関数として自然に定義される

10−5(11)1e-06
0.0001. 線型方程式の解法には許容誤差

0.001

0.01
frequency

0.1

1の Generalized minimal residual method (GMRES)(12)を用い, 右前処理の point Jacobi 法で前処理を行った. また
本稿の数値計算は, 導入した離散化法によって得られる数値

Fig. 3

Relative error for several frequencies withthe mesh of 1536 quadrilateral elements.

解の精度を確認することが主目的であるため, 簡単のために係数行列の計算の高速化アルゴリズムを用いなかった. すな

上記したとおり, 本稿の数値計算では係数行列計算の高速

わち線型方程式の係数行列は全ての要素を計算, 保存し, 係化を特に行っていないため, 全体の計算時間はほぼ係数行列数行列の乗算はこの保存した行列を用いて行った.の計算に要する時間に一致する. これより各数値解法に要す

Fig. 2, 3 は立方体形状の散乱体を, それぞれ 384 要素, 1536

る計算時間は周波数にほとんど依存しないと考えられる. 実

要素の四角形で分割したメッシュを用いて, 上記の三解法で際各数値解法の計算時間は入射波によらずほぼ一定の値を示した. Table 1 に周波数 ω = 0.1 のときの各数値解法の計算時た高速解法や Calderon の式に基づく前処理 (15) を導入する間を示す. Hdiv 内積を用いた解法では従来法と比べ, 要することで計算時間を短縮し, より複雑な問題に対して本稿で提Table 1

Computational time (s) with ω = 0.1

案した手法が有効であることを示すことが課題である. 特に高速解法を導入することで, 本稿で扱った数値例よりもさら

Number of elements

384

1536

DOF

1536

6144

Hdiv scalar product, c = c0424.690

6778.867

Hdiv scalar product, c = 1423.741

6789.911

Conventional method

401.8516412.467にメッシュの一辺の代表長 h の小さい問題を解くことができるようになると考えられる. これを用いて, 本稿では扱えなかった h の提案手法の精度に与える影響を調べることも課題の一つである.

参考文献(1) JD. Joannopoulos, RD. Meade, and JN. Winn. PhotonicCrystals. Princeton University Press, Princeton, NJ,

計算時間がやや多くなっているが, その増加量は自由度 1536

USA, 1995.のとき 5.68%, 自由度 6144 のとき 5.71%であり, Hdiv 内積を(2) JB. Pendry. Negative refraction makes a perfect lens.

用いた解法で新しく必要となる項の計算にかかる時間はそれ

Physical Review Letters, Vol. 85, pp. 3966–3969, 2000.

ほど多くないことがわかる.

(3) Weng Cho Chew, E Michielssen, JM Song, and JM Jin.

また Table 2 には周波数 ω = 0.1 のときの各数値解法における GMRES が要した反復回数を示す. 提案手法の反復回数

Table 2Iteration number with ω = 0.1

Fast and eﬃcient algorithms in computational electromagnetics. Artech House, Inc., 2001.
(4) Francesco Andriulli, Kristof Cools, Ignace Bogaert,and Eric Michielssen. On a well-conditioned electric

Number of elements

384

1536

DOF

1536

6144

tries. IEEE Transactions on Antennas and Propagation,

Hdiv scalar product, c = c041

81Vol. 61, No. 4, pp. 2077–2087, 2013.

Hdiv scalar product, c = 14079

Conventional method

66

124

field integral operator for multiply connected geome-(5) Zhi Guo Qian and Weng Cho Chew. An augmentedelectric field integral equation for high-speed interconnect analysis. Microwave and Optical Technology Letters, Vol. 50, No. 10, pp. 2658–2662, 2008.

は従来法と比較しやや少ない値を示しており, Hdiv (Γ) 内積を用いることで係数行列の条件が悪化することはないことが(6) Jianfang Zhu, Saad Omar, Wenwen Chai, and Dan
Jiao. A rigorous solution to the low-frequency breakdown in the electric field integral equation. In Antennas

わかる.

and Propagation (APSURSI), 2011 IEEE International
6. 結論

Symposium on, pp. 3214–3217. IEEE, 2011.

電磁波動散乱問題に対する PMCHWT 定式化による境界要素法において, Hdiv (Γ) 内積を用いた Galerkin 法によって特に周波数が小さい問題における精度を大きく改善できるこ(7) J.C. Nédélec.

Acoustic and Electromagnetic Equa-tions: Integral Representations for Harmonic Problems.
Springer Verlag, 2001.とを示した. また数値計算例を通して, 新しい解法では僅かな計算時間の増加で従来法の解の精度を大きく改善できるこ(8) Mei Song Tong, Weng Cho Chew, Barry J Rubin, Jason D Morsey, and Lijun Jiang.とを示した.
今後の課題として, 周波数 ω が比較的大きい値を持った問題に対して, 本稿で提案した Hdiv (Γ) 内積を用いた手法の適用可能性, 特に精度改善の効果があるかを調べることが挙げられ

On the dual ba-

sis for solving electromagnetic surface integral equations. IEEE Transactions on Antennas and Propagation, Vol. 57, No. 10, pp. 3136–3146, 2009.

る. 本稿では Hdiv (Γ) 内積中に現れる定数 c を ωh ≪ 1 の仮定(9) S. Rao, D. Wilton, and A. Glisson. Electromagneticの下決定し, ωh の値が小さい問題において提案手法が精度をscattering by surfaces of arbitrary shape. IEEE Trans-改善することのみを確認したが, 定数 c の値を適切に変更する

actions on Antennas and Propagation, Vol. 30, No. 3,

ことで様々な ω に対して効果的な数値解法を構成できること

pp. 409–418, 1982.

が期待される. また本稿は, 簡単な例を通して提案手法の効果(10) A. Buﬀa and S.H. Christiansen. A dual finite elementを確認するのみに留まった. 今後は Fast multipole method(13)complex on the barycentric refinement. Mathematics of(14)(FMM) や Adaptive cross approximation(ACA) といっ

Computation, Vol. 76, pp. 1743–1769, 2007.

(11) K. Niino and N. Nishimura. Calderón’s preconditioning
approaches for PMCHWT formulations for Maxwell’sequations.

International Journal of Numerical Mod-elling: Electronic Networks, Devices and Fields, Vol. 25,
pp. 558–572, 2012.
(12) Y. Saad.

Iterative Methods for Sparse Linear Sys-tems. Society for Industrial and Applied Mathematics,
Philadelphia, PA, USA, 2003.
(13) N Nishimura. Fast multipole accelerated boundary integral equation methods. Applied Mechanics Reviews,Vol. 55, p. 299, 2002.
(14) Sergej Rjasanow and Olaf Steinbach. The fast solution of boundary integral equations. Mathematical andanalytical techniques with applications to engineering.
Springer, New York, 2007.
(15) O. Steinbach and W.L. Wendland. The constructionof some eﬃcient preconditioners in the boundary element method. Advances in Computational Mathematics, Vol. 9, No. 1, pp. 191–216, 1998.