周波数を与えて伝搬定数を求める周期構造導波路のFEM表示式（２） －境界節点界のみを用いる方法－

周期構造導波路のFEM固有値解析で、界分布Φを直接解く方法でも周波数を与えて伝搬定数を求めることができましたが、この方法は問題を更に縮小できることが分かりました。

前回：

　　周波数を与えて伝搬定数を求める周期構造導波路のFEM表示式

　　　http://ryujimiya.hatenablog.com/entry/2013/01/30/014612

前回導出したFEM表示式をばらすとλ = exp( - jβd) として、

    [P00] {Φ}0 + [P01] {Φ}1 + λ [P02] {Φ}1 = {0}            (1)

    [P10] {Φ}0 + [P11] {Φ}1 + λ [P12] {Φ}1 + (1/λ) [P20] {Φ}0 + (1/λ) [P21] {Φ}1 + [P22] {Φ}1 = {0}          (2)

となりますが、ここで(1)より

    {Φ}0 = - [P00]-1 ( [P01] + λ [P02] ) {Φ}1                  (3)

(3)を(2)に代入して、{Φ}1だけの式にすると

    {  [C] + λ [M] + (1 / λ) [K]  }  {Φ}1 = {0}                    (4)    

           [C] =  - [P10] [P00]-1 [P01] + [P11] - [P20] [P00]-1 [P02] + [P22]

           [M] =  - [P10] [P00]-1 [P02] + [P12]

           [K] =  - [P20] [P00]-1 [P01] + [P21]

最終的な式は、λに関する２次の固有方程式になっているので、線形化して解くことができます。

この方法では固有方程式の次元が境界上の界のみなので、高速に解くことができます。

[P00]-1の逆行列の計算のコストが気になるところですが、私のPCではそれほどコストは高くなかったです。

それよりは、求めた固有ベクトル{Φ}1から内部の界ベクトル{Φ}0を求める処理の方が計算時間がかかりました。

これは、(3)にλが含まれるため固有値毎に毎回行列を作成する必要があるため、全固有値の界分布を求めようとすると時間がかかるわけです。

ただ、必要な固有値（つまり、伝搬モード）についてのみ全領域の固有ベクトルを計算するようにすると、それほど計算時間はかからないのでのでトータルの計算時間はこちらの方が領域全体の界に関する固有値問題を直接解く方法より短くなりました。

 

計算結果

　前回と同じエアホール型三角形格子フォトニック結晶欠陥導波路の伝搬定数を境界節点界のみを用いる方法で求めてみました。

 　エアホール半径r = 0.30a、基板の屈折率n = 2.76、60°三角形格子

前回の結果とよく一致しています。また、固有値問題を縮小したためか分かりませんが前回計算に失敗した周波数でもきちんと解が得られました。