光柵衍射問題數(shù)值計算的最小二乘方法

欒 天，王玉潔，鄭恩希，3

（1.北華大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院，吉林 吉林132033；2.吉林大學(xué) 數(shù)學(xué)研究所，長春130012；3.大連海事大學(xué) 理學(xué)院，遼寧 大連116026）

1 光柵衍射問題

在電磁場散射理論中，散射體是周期結(jié)構(gòu)表面的散射問題稱為光柵衍射［1］.光柵衍射在微光學(xué)領(lǐng)域（如光譜分析、納米尺度光學(xué)原件設(shè)計和光纖通信）應(yīng)用廣泛，關(guān)于這類問題的數(shù)值研究目前已取得了一些結(jié)果［2－4］，通常采用的數(shù)值方法為有限元方法［5－7］，適用于一般形狀的光柵結(jié)構(gòu).但當(dāng)問題的波數(shù)較大時，有限元方法將帶來較大的計算量［8］，不再適合實際應(yīng)用.為此，本文針對光柵衍射問題，提出一種更有效的算法——最小二乘方法，該方法不僅適用于一般形狀的衍射光柵，而且能克服大波數(shù)帶來的困難，應(yīng)用過程簡單，所需剖分單元少，收斂速度快。

考慮一維光柵的時諧衍射問題.設(shè)光柵的表面為Γ，周期為d.帶狀區(qū)域S定義為

設(shè)Γ0?S為Γ在帶狀區(qū)域S內(nèi)的一個周期，Γ0將S分割為兩個部分，上半部分記為E.入射平面波為uI＝exp｛iαx－iβy｝，其中α＝ksinθ，β＝kcosθ，k為波數(shù)，－π／2＜θ＜π／2為入射角.在區(qū)域E 中，全場u滿足Helmholtz方程：Δu＋k2u＝0.考慮滿足擬周期條件的解，即u滿足

這里u（x，y）exp｛－iαx｝在x方向是以d為周期的周期函數(shù).為了保證數(shù)學(xué)問題的解存在唯一，并且符合物理要求，還要求衍射場ud＝u－uI在無窮遠處滿足有界外行平面波條件，即衍射場在無窮遠處僅由有限多個向外傳播的平面波構(gòu)成.

一維光柵時諧衍射問題的數(shù)學(xué)模型為：給定平面入射波uI，求擬周期解u，滿足如下邊值問題：

且衍射場ud＝u－uI滿足有界外行平面波條件，這里a，b∈?不同時為零.當(dāng)a＝0，b≠0時，對應(yīng)Dirichlet邊界條件；當(dāng)a≠0，b＝0時，對應(yīng)Neumann邊界條件；當(dāng)a≠0，b≠0時，對應(yīng)Robin邊界條件.

2 最小二乘方法

將一個周期內(nèi)的計算區(qū)域進行簡單剖分，在每個有界子域中選用平面波函數(shù)近似解u的局部性態(tài).關(guān)于解在無窮遠處的性態(tài)，使用Rayleigh展開的前有限項截斷去近似，這種近似可以自然滿足有界外行平面波條件和擬周期條件.

考慮 子 區(qū) 域 Ej（j＝1，2，…，s）.由 于Helmholtz方程的解u∈H1（Ej）可以利用平面波函數(shù)去逼近［9］，因此可采用平面波函數(shù)近似解u的局部性態(tài).

在每個Ej內(nèi)部選取點xj＝（xj，yj）（j＝1，2，…，s），并選取Nj個方向d（j）l（l＝1，2，…，Nj），定義局部近似空間

圖1 符號示意圖Fig.1 Diagram of notations

在子區(qū)域Es＋1中，u的 Rayleigh展式［2］為

其中：

因此，可以采用Rayleigh展開的有限項截斷作為解u的近似.在Es＋1中，定義近似空間

結(jié)合上述所有近似空間Vj（j＝1，2，…，s＋1），定義試探函數(shù)空間V為

于是，可定義誤差匹配泛函為

其中［v］表示函數(shù)v在子區(qū)域相交邊界處的躍度，定義為

利用對偶技巧可以得到最小二乘方法的一個基本估計：

命題1 若k2≠（αn＋α）2，n∈?，則存在一個與u和Nj無關(guān)的常數(shù)C，使得

命題1的證明過程與文獻［10］中定理3.1和文獻［11］中定理2.1的證明完全類似，故略.命題1表明，對任意非共振波數(shù)k，J（uN）1／2控制著解的內(nèi)部誤差，因此可以用于判斷算法的收斂性.

3 數(shù)值模擬

下面通過數(shù)值模擬驗證算法在計算光柵衍射問題時的有效性.數(shù)值實驗均使用Matlab軟件實現(xiàn).

為簡單，取Nj＝p（j＝1，2，…，s＋1），p∈?，即每個單元上選取相同數(shù)目的平面波函數(shù)，且平面波的方向按如下方式選取：d（j）l＝（cosθl，sinθl），θl＝2π（l－1）／p，l＝1，2，…，p.

例1 直線光柵Γ＝｛（x，y）∈?2，y＝0｝，入射波為平面波uI＝exp｛iαx－iβy｝.取波數(shù)k＝50，光柵周期d＝2，入射角θ＝π／4，a＝0，即Dirichlet邊界條件.將計算區(qū)域｛（x，y）；－1＜x＜1，y＞0｝剖分，選取點x1＝（－0.5，0.25），x2＝（0.5，0.25），如圖2所示.

由于例1中問題存在真解u＝uI－exp｛iαx＋iβy｝，因此本文計算了真解和數(shù)值解在區(qū)間Ω＝（－1，1）×（0，1）上的L2誤差，并分析了L2誤差和泛函J（uN）1／2的收斂性.數(shù)值結(jié)果表明，誤差隨基底數(shù)目的增加而快速減少，而且L2誤差可以被J（uN）1／2所控制，如圖3所示.

圖2 直線光柵與區(qū)域剖分Fig.2 Straight line grating and domain decomposition

圖3 L2 誤差與J（uN）1／2收斂結(jié)果Fig.3 L2 error and convergence result of J（uN ）1／2

圖5為泛函J（uN）1／2的收斂性結(jié)果.由圖5可見，隨基底數(shù)目的增加，泛函值隨之快速衰減.當(dāng)基底數(shù)目達到收斂性要求時，收斂速度也很快.此外，由于L2誤差可以被泛函J（uN）1／2所控制，所以即使在沒有真解的情況下，仍可以判定算法是關(guān)于基底數(shù)目p收斂的.

圖4 曲線光柵與區(qū)域剖分Fig.4 Curve grating and domain decomposition

圖5 J（uN）1／2收斂結(jié)果Fig.5 Convergence result of J（uN ）1／2

由例1和例2可見，本文算法在處理光柵衍射問題時是高效的.一方面，算法僅需較少的剖分單元，從而減小了計算量；另一方面，當(dāng)波數(shù)較大時，算法同樣可以達到較好的精度.

［1］MENG Pin－chao，JIANG Zhi－xia， LI Yan－zhong.Electromagnetic Scattering of Diffractive Grating in Homogeneous Medium ［J］.Journal of Jilin University：Science Edition，2012，50（6）：1151－1155.（孟品超，姜志俠，李延忠.均勻介質(zhì)中衍射光柵的電磁散射 ［J］.吉林大學(xué)學(xué)報：理學(xué)版，2012，50（6）：1151－1155.）

［2］BAO Gang，Cowsar L，Masters W.Mathematical Modeling in Optical Science［M］.Philadelphia：Frontiers in Applied Mathematics，2001.

［3］YIN Wei－shi，ZHANG De－yue，MA Fu－ming.Numerical Calculation of the Scattering Problem for Grating by Integral Equation Method［J］.Journal of Jilin University：Science Edition，2009，47（6）：1112－1121.（尹偉石，張德悅，馬富明.光柵散射問題數(shù)值計算的積分方程方法 ［J］.吉林大學(xué)學(xué)報：理學(xué)版，2009，47（6）：1112－1121.）

［4］LUAN Tian，MA Fu－ming.Well－Posedness of Anisotropic Layers Scattering above Rough Surfaces［J］.Journal of Jilin University：Science Edition，2012，50（2）：213－218.（欒天，馬富明.粗糙曲面上各向異性介質(zhì)層散射問題的適定性 ［J］.吉林大學(xué)學(xué)報：理學(xué)版，2012，50（2）：213－218.）

［5］BAO Gang.Numerical Analysis of Diffraction by Periodic Structures：TM Polarization［J］.Numer Math，1996，75（1）：1－16.

［6］BAO Gang，CAO Yan－zhao，YANG Hong－tao.Numerical Solution of Diffraction Problems by a Least－Squares Finite Element Method［J］.Math Methods Appl Sci，2000，23（12）：1073－1092.

［7］BAO Gang，CHEN Zhi－ming，WU Hai－jun.Adaptive Finite－Element Method for Diffraction Gratings［J］.J Opt Soc Am A，2005，22（6）：1106－1114.

［8］Deraemaeker A，Babu?ka I，Bouillard P.Dispersion and Pollution of the FEM Solution for the Helmholtz Equation in One，Two and Three Dimensions［J］.Int J Numer Meth Eng，1999，46（4）：471－499.

［9］Moiola A，Hiptmair R，Perugia I.Plane Wave Approximation of Homogeneous Helmholtz Solutions ［J］.Z Angew Math Phys，2011，62（5）：809－837.

［10］ZHENG En－xi.Application of Least－Squares Non－polynomial Finite Element Methods in Several Scattering Problems［D］.Changchun：Jilin University，2012.（鄭恩希.非多項式最小二乘有限元法在幾種散射問題中的應(yīng)用 ［D］.長春：吉林大學(xué)，2012.）

［11］ZHENG En－xi，MA Fu－ming，ZHANG De－yue.A Least－Squares Non－polynomial Finite Element Method for Solving the Polygonal－Line Grating Problem ［J］.Journal of Mathematical Analysis and Applications，2013，397（2）：550－560.