SENAC Vol. 55, No. 1（2022. 1）― 3 ―[共同研究成果]［共同研究成果］

ナノアンテナの数値解析のためのブロックモーメント法の高速化
今野 佳祐 † , 陳 強 †東北大学大学院工学研究科 通信工学専攻 †まえがき1

近年，筆者らはナノアンテナの数値解析のためのブロックモーメント法の高速化を行ってきた．昨年度の報告では，ブロックモーメント法を用いてインピーダンス行列を数値計算する部分のベクトル化を実施し，その一部を高速化することに成功した [1]．その一方で，インピーダンス行列の数値計算を行う部分のうち，指数積分の数値計算のベクトル化が困難であったため，インピーダンス行列の数値計算アルゴリズム全体のベクトル化には成功していなかった．指数積分の数値計算では，精度のよい結果が得られるアルゴリズムが引数の値や積分路によって異なる．そのため，引数の値や積分路によっては必ずしもベクトル化しやすい数値計算アルゴリズムを選択できないことがあり，ベクトル化には工夫が必要である．

本報告では，昨年度に引き続きブロックモーメント法のベクトル化を行い，インピーダンス行列の数値計算アルゴリズム全体のベクトル化を行った結果を示す．ベクトル化が困
難な指数積分の数値計算の部分は，精度を損なわない範囲でベクトル化に有利なアルゴリズムに変更し，数値計算の高速化を図った．数値計算例を示し，高速化アルゴリズムの有効性を明らかにする．

体積積分方程式を用いたブロックモーメント法

22.1

線状モノポールセグメント間の自己・相互インピーダンス

さて, 図 1 の観測セグメント上の電流分布を表す基底関数として, 以下の区分正弦関数を利用する [2]，[3]．

sinh jk0 (ea − ea1 ).
sinh jk0 |ea2 − ea1 |(1)sinh jk0 (sa2 − sa ).
fs (s ) =sinh jk0 |sa2 − sa1 |(2)

fe (ea ) =同様に，波源セグメント上の電流分布を表す基底関数として，以下の区分正弦関数を用いる．

a― 4 ―

SENAC Vol. 55, No. （1 2022. 1）

次に, (1) 式を連続の式に代入することによって, 観測セグメント上の線電荷密度を以下のように得る.

∇ · fe (ea ) + jωρe (ea ) = δ(ea − ea2 ){}1sinh jk0 (ea − ea1 )
aaa−jk0+ δ(e − e2 ) .
ρe (e ) =jω
sinh jk0 |ea2 − ea1 |(3)同様に, (2) 式を連続の式に代入することにより, 波源セグメント上の線電荷密度を以下のように得る.

∇ · fs (sa ) + jωρs (sa ) = −δ(sa − sa1 ){}1sinh jk0 (sa2 − sa )aaa
jk0− δ(s − s1 ) .
ρs (s ) =jω
sinh jk0 |sa2 − sa1 |(4)ここで, 以下に示すセグメント間の自己・相互インピーダンスの一般的な表示式に (1)-(4)

式を代入して整理すると, 端電荷を持つ線状モノポールセグメント間の自己・相互インピー

ダンス表示式が, 単積分の形で (5) 式のように得られる [4].
∫ ea2 ∫ sa2∫ a∫ a−jk0 Rijjωµ0 e2 s2jω
e−jk0 Rij a aaa eaa
Zij =fe (e ) · fs (e )ds de +ρe (ea )ρs (sa )ds de
4π ea1 sa1Rij
4πε0 ea1 sa1Rij= C A I A + CB I B + CC I C − C D I D .

(5)√ただし, (5) 式中の CA ∼ CD , IA ∼ ID , そして Rij はそれぞれ，Rij = (ea − sa )2 + (sb2 ) + (sc2 )2 ,22 ∑∑aae−jk0 (ne1 +ms2 ) mnIpq , CB = j sin k 15sa −sa , IB =CA = sin k ea −ea7.5sin k sa −sa , IA =0| 20| 20| 21|1|1|
p=1 q=122∑∑
aa ∫ v(sa ,ea ) −z
IBp , IBp = −mejk0 (me2 −ms2 ) v(sa2,ea2) e z dz, CC = j sin k 15ea −ea , IC =ICq , ICq =1 20| 21|p=1
q=1−jk0 R21aa ∫ w(sa ,ea ) −z
30, ID = e R21である．なお，IA の表示式に含−nejk0 (−ne1 +ns1 ) w(sa1,ea2) e z dz, CD = jk0
まれる Ipq は以下のようになる．
∫ a∫ aaa1 + mn e2 s2 ejk0 (ms +ne −Rij ) a a
Ipq = jk0
ds demnRijea
sa11∫ v(sa1 ,eah ) −z
∫ w(sah ,ea2 ) −z2∑eea (m+n)
hjk0 ea(n+m)
jks0dz − e hdz=(−1) {−e ha ,ea )a ,ea )zzv(sw(s1 hh 1h=1
∫ x(sah ,ea2 ) −z
∫ y(sah ,ea2 ) −z
ee+ ejk0 ju0dz + e−jk0 ju0dz}.
aazz
x(say(sah ,e1 )h ,e1 )11(6)ただし，v(sai , eah ) = jk0 (−msai +meah +Rhi ), w(sah , eai ) = jk0 (−neai +msah +Rih ), x(sah , eai ) =√, d=
jk0 (−msah −neai +Rih +ju0 ), y(sah , eai ) = jk0 (−msah −neai +Rih −ju0 ), u0 = d 1+mn
1−mn√
(sb2 )2 + (sc2 )2 , m = (−1)p−1 , n = (−1)q−1 であり，(6) 式において, 1 + mn = 0 ならIpq = 0, 1 − mn = 0 なら最後の 2 項を無視する. Ipq ，IB ，IC の表示式において, z を含む全ての積分は 2 つの指数積分の差を取ることで数値的に積分できる．なお，モノポールセグメントの配置が垂直の場合も同様の定式化が必要だが, そちらは省略する．

Ipq ，IB ，IC の表示式において，指数積分の数値計算が最も時間のかかる部分である．

ナノアンテナの数値解析のためのブロックモーメント法の高速化

― 5 ―

eai

d

( s2b )2  ( s2c ) 2 (e2a , 0,0)bbcc
ass,ss(e1 , 0, 0)  12 12ebOdj( s2a , s2b , s2c )( s1a , s1b , s1c )ec
図 1: 端電荷を持つ一対の線状モノポールセグメントの一例と座標系.

2.2

ブロックモノポールセグメント間の自己・相互インピーダンス
2 weceale iOec( s2b )2  ( s2c ) 2 (e2a , 0,0)  dbcc
s1b s,ss(e1a , 0, 0) 
2 12

2 webebc
2ws( s2a , s2b , s2c )dls j( s1a , s1b , s1c )2 wsb

図 2: 端電荷を持つ一対のブロックモノポールセグメントの一例と座標系.
図 2 に, 端電荷を持つ一対の平行なブロックモノポールセグメントを示す. 観測・波源

モノポールセグメントにおける電流の向き, 端電荷の符号と位置は前項の図 1 と同じである. 観測ブロックモノポールセグメントの電流分布を表す基底関数は, 以下のように定義される.

fe (ea ) =1 sinh jk0 (ea − ea1 ).
4web wec sinh jk0 |ea2 − ea1 |(7)

1 sinh jk0 (sa2 − sa ).
4wsb wsc sinh jk0 |sa2 − sa1 |(8)

同様に, 波源モノポール上の電流分布を表す基底関数は以下の式で定義される.

fs (sa ) =これらの基底関数は, 先頭の係数部分を除けばそれぞれ (1) 式, (2) 式と一致する. 従って,ブロックセグメント中央を通る線状モノポールセグメント間の相互インピーダンス Zij を(5) 式によって求め, Zij を 4 重積分すると, ブロックモノポールセグメント間の自己・相互インピーダンス表示式が以下のように得られる.
∫ web ∫ wec ∫ sb1 +wsb ∫ sc1 +wsc
1blockZij dsc dsb dec deb .
Zij =16web wec wsb wsc −web −wec sb1 −wsb sc1 −wsc(9)

Zij の計算には単積分が含まれるため, (9) 式は 5 重積分を含む表示式である. この 5 重積

分は，座標変換を用いて 3 重積分の和 (互いに垂直なブロックセグメントであれば 4 重積

― 6 ―

SENAC Vol. 55, No. （1 2022. 1）分) に変換することができ，数値計算の収束性および計算時間が改善する [5]．詳細は省略

するが，本報告ではこの 3 重積分の表現を用いて数値計算を行う．ナノアンテナの数値解
析では，これらの式を用いて誘電体ブロックセグメント間の自己・相互インピーダンス行列要素を計算し，得られた行列方程式を解くことで電流を求める．

2.3

指数積分のベクトル化

誘電体ブロックセグメント間の自己・相互インピーダンスの数値計算において，指数積分の数値計算には最も長い計算時間がかかる．したがって，指数積分の数値計算のベクトル化は自己・相互インピーダンスの数値計算の高速化を図る上で必須であるが，NECSX-Aurora TSUBASA 向け数学ライブラリ（NEC Numeric Library Collection）には複素数を引数とする指数積分のライブラリがない [6]．そこで，本報告では指数積分の数値計算のベクトル化を独自に試みた．本報告での指数積分の数値計算に用いたサブルーチンは，[7] からダウンロードしたものであるので，そのベクトル化を実施した．

指数積分はいわゆる複素積分である．被積分関数は原点に極を持ち，積分路がブランチカット (=負の実軸) を横切ると異なるリーマン面へと移動する多価関数である．したがって，積分路が極の近くを通ったり，ブランチカットを横断するときの積分には注意が必要である．このような指数積分を数値的に求める方法には，漸近展開や級数展開，数値積分の他，連分数展開を用いる方法が挙げられる．これらの中で，漸近展開や級数展開，数値積分を用いたアルゴリズムのベクトル化は容易であるが，連分数展開を用いたアルゴリ

ズムはループ内の変数に依存関係ができてしまうのでベクトル化が困難である．したがって，指数積分の数値積分の高速化だけを考えるのであれば，連分数展開以外のアルゴリズムを用いればよいことになる．一方で，数値計算の精度を考えると，引数の値によって数値計算のアルゴリズムを変えるのが最もよいことになり，引数の値によっては連分数展開を用いることは避けられない．

本報告では，以上の議論を踏まえ，引数の値の範囲と取り得る積分路によって指数積分の数値計算のアルゴリズムを変える．互いに平行な誘電体ブロックセグメント間の自己・相互インピーダンスの数値計算では，指数積分の引数は純虚数になり，積分路が原点を横断することがない．そこで，このときの指数積分の数値計算は全て数値積分で実施する．
一方，互いに垂直な誘電体ブロックセグメントの自己・相互インピーダンスの数値計算

では，指数積分の積分路が原点あるいはブランチカットを横断することがある．そこで，このときは場合分けを行い，その結果に基づいて指数積分の数値計算アルゴリズムを変える．指数積分の引数が複素数のとき，被積分関数の特異性が強い可能性がある (要するに，積分路が極やブランチカット付近を通る可能性がある)．したがって，このときは精度を優先し，数値計算は連分数展開と級数展開を用いて行う．指数積分の引数が純虚数のとき，被積分関数の特異性が弱ければ (要するに，積分路が極を横断していなければ) 高ナノアンテナの数値解析のためのブロックモーメント法の高速化― 7 ―速化を優先し，その数値計算は数値積分で行う．反対に，被積分関数の特異性が強ければ精度を優先し，その数値計算は連分数展開と級数展開で行う．このように，引数の値の範囲と取り得る積分路によって指数積分の数値計算のアルゴリズムを変えることで，精度を損なうことなくベクトル演算率の向上と高速化を図る．なお，数値計算に級数展開や連分数展開を用いる箇所にもループアンローリングを施し，可能な限りベクトル演算率を向上させた．

3数値計算結果
100000

CPU time [sec.]

10000
1000
1000.001N2101110100
1000
Number of unknowns N10000

図 3: 数値計算時間．

銀でできた 100 × 100 × 100 nm3 の立方体を垂直入射の TM 波で励振し，その電流を数値的に求めた．銀の複素比誘電率は Drude モデルによって与え，周波数は 100 THz とした [8]．ナノダイポールアンテナはブロックダイポールおよびブロックモノポールセグメントに分割した．

このときのセグメント数 N に対するトータルの数値計算時間を図 3 に示す．数値計算時間は概ね N 2 に比例しており，これはインピーダンス行列の数値計算時間のオーダーと一致する．一般的に，モーメント法の数値解析では行列方程式を解くための数値計算時間が N 3 に比例して増加するので，トータルの数値計算時間は N 3 に比例して増加することが多い．一方，今回の数値計算では，インピーダンス行列の数値計算に多重積分が含まれており，そのせいでインピーダンス行列の数値計算がトータルの計算時間の大部分を占めている．その結果，トータルの数値計算時間が N 2 に比例して増加している．このよう― 8 ―

SENAC Vol. 55, No. （1 2022. 1）に，誘電体の数値解析では，自己・相互インピーダンスの数値計算の高速化が必須であることが分かる．

プログラム全体のベクトル演算率は 95.3%であり，指数積分の箇所以外のベクトル化を実施した昨年度のプログラムのベクトル演算率 2.7%と比較して大幅にベクトル演算率が向上した．自己・相互インピーダンス 1 つあたりの数値計算時間は 0.0016 sec. であり，昨年度のプログラムの数値計算結果の 0.43 sec. と比較して 200 倍以上の高速化を実現した．

4

まとめ

本報告では，ナノアンテナの数値解析のための体積積分方程式を用いたモーメント法の高速化を行い，指数積分の数値計算の箇所のベクトル化に成功した．その結果，プログラム全体の高速化が図れ，数値計算の大幅な高速化に成功した．今後は，様々なモデルで本プログラムの数値計算精度を検討し，アンテナの設計にも応用していく予定である．

謝辞

本研究成果は, 東北大学サイバーサイエンスセンターのサポートを受けたものである．

スタッフをはじめ，関係各位に深謝する．

参考文献[1] 今野 佳祐，陳 強，“ナノアンテナの数値解析のためのブロックモーメント法の高速化，” 東北大学大規模科学計算システム広報 SENAC vol. 54, no. 1, pp.25-31, Jan.

2021.

[2] J.H. Richmond and N.H. Geary, “Mutual impedance of nonplanar-skew sinusoidaldipoles,” IEEE Trans. Antennas Propag., vol.AP-23, no.3, pp.412-414, May 1975.
[3] 陳 強，ザイ フイチン，袁 巧微，澤谷 邦男，“誘電体に対するガラーキンモー
メント法 -端部電荷を考慮した直方体モノポール間の自己・相互インピーダンスの単積分化-，” 信学論 (B)，Vol. J91-B，No. 9，pp. 926-939，2008 年 9 月．

[4] M.A. Tilston and K.G. Balmain, “On the suppression of asymmetric artifacts arisingin an implementation of the thin-wire method of moments,” IEEE Trans. Antennas
Propag., vol.38, no.2, pp.281-285, Feb. 1990.
[5] A. Köksal and J.F. Kauﬀman, “Mutual impedance of parallel and perpendicularcoplanar surface monopoles,” IEEE Trans. Antennas Propag., vol.39, no.8, pp.12511256, Aug. 1991.

ナノアンテナの数値解析のためのブロックモーメント法の高速化― 9 ―[6] NEC Numeric Library Collection 2.3.0 ユーザーズガイド，https://www.hpc.nec/

documents/sdk/SDK_NLC/UsersGuide/asl/c/ja/index.html, 2021 年 8 月 16 日ア

クセス．

[7] Special Functions, “https://people.sc.fsu.edu/~jburkardt/f77_src/special_functions/special_functions.html,” 2021 年 8 月 6 日アクセス．

[8] M. Nafari, and J. M. Jornet, “Modeling and performance analysis of metallic plasmonic nano-antennas for wireless optical communication in nanonetworks,” IEEE
Access, vol. 5, pp. 6389-6398, 2017.