趙艷敏,何滄平

矢量波動方程的時空高階方法

趙艷敏1,何滄平2,3

(1.許昌學院數(shù)學科學學院,河南許昌461000; 2.中國科學院數(shù)學與系統(tǒng)科學研究院計算數(shù)學與科學工程計算研究所,北京100190;3.中國科學院研究生院,北京100190)

文章將Gauss-Lobatto-Legendre多項式的高階矢量譜元方法應用于矢量波動方程.由于矢量波動方程可以表示為一個無窮維Hamilton系統(tǒng)且經空間上的有限元方法離散后是一有限維Hamilton系統(tǒng),利用4階辛分塊的Runge-Kutta方法來求解該有限維Hamilton系統(tǒng),以期保持系統(tǒng)整體的能量和結構.

矢量波動方程;矢量譜元方法;Hamilton系統(tǒng);辛分塊的Runge-Kutta方法

0 引言

由于電磁場計算對生物電磁學、醫(yī)學、微波遙感、通訊技術和納米電子學等眾多領域的重要影響,對電磁場問題如何設計數(shù)值算法成了人們越來越關注的課題.矢量波動方程就是其中重要的一類.求解此類方程,工程界用得比較多的方法是有限差分方法[1-3].隨后,有限元方法憑借其區(qū)域靈活性優(yōu)勢也越來越多地被應用于計算電磁場問題[4-7].最近,柳清伙等學者又提出了利用譜元方法與傳統(tǒng)的Runge-Kutta方法相結合來求解時域Maxwell方程[5].矢量波動方程可以表示成一個無窮維Hamilton系統(tǒng),因此它的解可以看做是泛函空間的Hamilton流,它在時間方向上是保整體能量和辛結構的.辛方法是唯一能精確保持連續(xù)Hamilton流內在性質的數(shù)值方法.基于此優(yōu)勢,辛算法是求解由無窮維Hamilton系統(tǒng)通過離散方法降維而成的有限維Hamilton系統(tǒng)的首選[8-13].本文對于矢量波動方程,給出了其在空間上經基于Gauss-Lobatto-Legendre多項式的高階譜元方法離散后的有限維Hamilton系統(tǒng);并給出了在時間方向上全離散求解該有限維系統(tǒng)的高階辛方法-辛分塊的Runge-Kutta方法.本文的內容安排如下:第一節(jié),給出了矢量波動方程及其無窮維Hamilton形式,并對空間半離散后的系統(tǒng)進行分析;第二節(jié),給出了空間離散所采用的高階矢量譜元方法及剛度矩陣、質量矩陣的一些計算技巧;第三節(jié)是對本文內容的總結.

1 矢量波動方程和高階辛算法

在各向同性、非耗散的線性介質中,考慮如下矢量波動方程:

其中ε和μ分別是介質的介電常數(shù)和磁導率.

設p為標量,q=(q1,q2)為向量,則有

引入變量F,(1)可以改寫為:

為簡單起見,設ε,μ為常數(shù),我們考慮的區(qū)域Ω中有有限種介質.上述矢量波動方程可以表示為如下無窮維的Hamilton系統(tǒng):

在空間方向,(2)式的弱形式為:找E,F∈H0(curl),?W∈H0(curl)使得

其中〈U,W〉=∫ΩU·WdΩ;H0(curl)={v∈(L2(Ω))2;curlv∈L2(Ω),n×v|?Ω=0}.令Γh為區(qū)域Ω的矩形正則剖分.設^K是ε-η平面的參考元,對剖分Γh而言,對于任意的K∈Γh,存在一個可逆映射FK:^K→K.

令Φi,i=1,2,…,Ne為^K上的基函數(shù).相應的有限元空間為(在第二節(jié)中,我們給出具體的有限元空間):

相應的離散問題是:對于(4),找E,F∈Vh,對任意的ˉΦj∈Vh,使得

顯然,在K上,E,F可以由基函數(shù)線性表出:

將它們代入(5)式,可得其矩陣形式:

總體質量矩陣M和總體剛度矩陣S分別由單元質量矩陣M(e)和單元剛度矩陣S(e)按照總體自由度編碼構成,其中M(e)和S(e)的元素為:

其中

由于原方程經有限元離散后為一有限維Hamilton系統(tǒng),我們用4階顯式辛分塊的Runge-Kutta (SPRK)方法對式(7)進行全離散求解.s階辛分塊的Runge-Kutta(SPRK)方法的系數(shù)可由圖1給出.

圖1 S階辛塊Runge-Kutta方法的Butcher表Fig.1 Table of Butcher for SPRK method

通常,它應用于常微分方程

結果為:

PRK方法的辛條件是:

(7)式是一個可分系統(tǒng),故PRK方法的辛條件僅有(8)式.

特別地,4階顯式分塊Runge-Kutta(SPRK)方法的Butcher表見圖2,其中γ1=

圖2 4階辛分塊Runge-Kutta方法的Butcher表Fig.2 Table of Butcher for 4th-order symplectic SPRK method

2 高階矢量譜元方法

我們采用基于N階Gauss-Lobatto-Legendre(GLL)多項式的矢量譜元來離散矢量波動方程.在第一節(jié)中可以看到,最后用辛方法求解方程(7)時,涉及的只是通過有限元離散后得到的總體質量矩陣和總體剛度矩陣.在本節(jié)中,我們介紹矢量的GLL譜元方法并在保證逼近精度的前提下,利用一些技巧在求總體質量矩陣和總體剛度矩陣及其逆矩陣時,盡可能地減少計算量.

N階Legendre多項式是

GLL點{ξj,j=0,1,…,N}是多項式(1-ξ2)L′N(ξ)的零點.

在參考元[-1,1]上,N階GLL基函數(shù)為:

顯然有φj(ξk)=δjk,?j,k=0,…,N.

在(7)式中,我們令Φi=Φξi+Φηi,其中Φξi=Φξrs=ξ^ φr(ξ)φs(η);Φηi=Φηrs=η^ φr(ξ)φs(η).

下面給出S(e),M(e)的具體表達.為了計算方便,我們將S(e)進行分塊處理:S(e)=[A B;C D]

各塊的元素分別為:

計算過程中,我們將用到GLL多項式φj在第i個GLL點的導數(shù)值,它即是矩陣～D的元素ˉDij:

令n=N+1,在參考元[-1,1]×[-1,1]上,我們易得:

為了節(jié)省計算時間、減少計算量,我們將采用一些技巧使得S(e),M(e)對角或塊對角.我們以矩陣A和M(e)的計算為例來展示:

單元質量矩陣M(e)的元素為

設K是矩形單元,則Jacobi陣G是一對角陣,即G=diag(g11,g22).

根據(jù)我們的自由度編碼,以i,j∈[1,n2]為例,即Φi=^ξ φr(ξ)φs(η),Φj=^ξ φr′(ξ)φs′(η),處理M(e)的過程如下:

S(e),M(e)已求出,則按照總體自由度編碼以及總體自由度與單元自由度的對照關系,容易得到總體剛度矩陣S及其逆陣S-1和質量陣M及其逆陣M-1.最后,利用第一節(jié)中所述的4階辛分塊的Runge-Kutta方法即得方程(7)全離散解.

3 總結

本文利用基于GLL多項式的矢量譜元空間離散方法結合時間方向上的辛離散方法來求解矢量波動方程.充分利用了譜元方法的高精度特點,并采用了一些數(shù)值處理技巧降低計算量從而解決了高階方法計算量大的問題;注意到矢量波動方程譜元半離散后的系統(tǒng)的結構特點,文中采用的高階辛格式在保證解精度的前提下,保持了系統(tǒng)整體的能量和結構.

[1] YEE K S.Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media[J]. I EEE Trans A nte Prop,1966,14:302-307.

[2] TAFLOVE A,HAGNESS S C.Computational Electrodynamics:The Finite-Difference Time-Domain Method3rd[M].Boston:A rtech House,2005,Boston.

[3] SHA W,HUANG Z X,WU X L,CHEN M S.Application of the Symplectic Finite-Difference Time-Domain Scheme to Electromagneticsimulation[J].J Comput Phys,2007,225:33-50.

[4] MONK P B.Analysis of A Finite Element Method for Maxwell’s Equations[J].S IA M J N umer A nal,1992,29:714-729.

[5] LIU Y X,L EEJ H,XIAO T,LIU Q H.A Spectral-Element Time-Domain Solution of Maxwell’s Equations[J].MicroOpti Technol Lett,2006,48:673-680.

[6] RIEBEN R,WHITE D,RODRIGUE G.High-Order Symplectic Integration Methods for Finite Element Solutions to Time Dependent Maxwell’s Equations[J].I EEE Trans A nte Prop,2004,52:2190-2198.

[7] J IAO D,J IN J M.Three-Dimensional Orthogonal Vector Basis Functions for Time-Domain Finite Element Solution of Vector Wave Equations[J].I EEE Trans on A nte and Prop,2003,51:59-66.

[8] FENG K.Collected Works of Feng Kang(II)[M].Beijing:National Defense Industry Press,1995.

[9] QIN M Z,ZHU W J.Construction of Higher Order Symplectic Schemes by Composition[J].Computing,1992,47:309-321.

[10] LU X W,SCHMID R.A Symplectic Algorithm for Wave Equation[J].Math Comput Simul,1997,43:29-38.

[11] TANG Y F,CAO J W,LIU X T,SUN Y C.Symplectic Methods for the Ablowitz-Ladik Discrete Nonlinear Schr?dinger Equation[J].J Phys A:Math Theor,2007,40:2425-2437.

[12] SUN G.Construction of High Order Symplectic Runge-Kutta Methods[J].J Comput Math,1993,11:250-260.

[13] ZHAO Y M,DAI G D,TANG Y F,LIU Q H.Symplectic Discretization for Spectral Element Solution of Maxwell’s E-quations[J].J Phys A:Math Theor,2009,42:325203.

Time-Domain High-Order Method for Vector Wave Equation

ZHAO Yan-min1,HE Cang-ping2,3
(1.School of Mathematical Sciences,Xuchang University,Xuchang461000,China; 2.L S EC,ICMS EC,Academy of Mathematics&S ystems Science,Chinese Academy of Sciences,Beijing100190,China; 3.Graduate School of Chinese Academy of Sciences,Beijing100190,China)

The high-order vector SEM based on Gauss-Lobatto-Legendre(GLL)polynomials is applied to vector wave equation.Since the vector wave equation can be denoted as an infinite-dimensional Hamiltonian system which will become a finite-dimensional Hamiltonian system by the finite element discretization on spatial direction,we adopt 4th-order symplectic partitioned Runge-Kutta(SPRK)method to solve the finitedimensional Hamiltonian system for conserving total energy and structure of the system.

vector wave equation;vector spectral element method;Hamiltonian system;symplectic partitioned Runge-Kutta method

O241

A

0253-2395(2010)03-0329-06

2010-01-11

國家自然科學基金(10672143;10872037);河南省教育廳自然科學基金(2009A110017)

趙艷敏(1979-),女,講師,博士,研究領域:辛算法、有限元理論及其應用.E-mail:zhaoym@lsec.cc.ac.cn