劉雙兵 陳海波 楊漢生

（巢湖學(xué)院物理與電子科學(xué)系，安徽 巢湖 238000）

Hamilton系統(tǒng)中Maxwell方程組的數(shù)值求解

劉雙兵 陳海波 楊漢生

（巢湖學(xué)院物理與電子科學(xué)系，安徽 巢湖 238000）

在詳細(xì)闡述Hamilton系統(tǒng)中的辛算法的基礎(chǔ)上，給出了Maxwell方程組的Hamilton的函數(shù)形式，將辛算法保持能量守恒和辛對稱性的思想應(yīng)用于Maxwell方程組的數(shù)值求解，結(jié)合傳統(tǒng)的時域有限差分（FDTD）法，得到了電磁場時間和空間的離散迭代公式，即辛?xí)r域有限差分（SFDTD）法。最后扼要地分析了該數(shù)值計算方法的穩(wěn)定性及數(shù)值色散性。

哈密爾頓系統(tǒng)；辛算法；時域有限差分；數(shù)值色散性

1 引言

近年來，隨著計算機(jī)性能的飛速發(fā)展和計算數(shù)學(xué)、計算物理中各種新型算法的出現(xiàn)，計算電磁學(xué)呈現(xiàn)出空前繁榮的局面。各種電磁場數(shù)值方法層出不窮，但這些方法面臨計算時間、存儲空間及計算精度等方面的困難，而且隨著人們對問題的物理本質(zhì)認(rèn)識的深入，意識到在追求算法高精度的同時，還應(yīng)力求保持原系統(tǒng)的內(nèi)在性質(zhì)。由于電磁場方程可以轉(zhuǎn)化為一無窮維Hamilton系統(tǒng)，而Hamilton系統(tǒng)具有一系列的內(nèi)在性質(zhì)，因而在對Hamilton系統(tǒng)的數(shù)值求解時，保持其內(nèi)在性質(zhì)就顯得尤為重要。辛算法正是保持Hamilton系統(tǒng)內(nèi)在性質(zhì)的一種新型數(shù)值方法，該算法在長時間的數(shù)值計算中，具有常見數(shù)值方法無可比擬的計算優(yōu)勢。

時域有限差分(FDTD)[1]法以其簡單直觀的特點已成為電磁學(xué)數(shù)值計算的一種常用方法，它直接求解依賴時間的Maxwell旋度方程，利用二階精度的中心差分近似旋度方程中的時間及空間微分算符，極易處理非均勻媒質(zhì)的情形。然而，它的計算精度相對較低，計算的時間步長與空間離散網(wǎng)格的大小必須滿足Courant-Fredrichs-Lewy(CFL)穩(wěn)定性條件，并且隨著計算時間步的增加，誤差將會累積。為了克服FDTD方法的這些缺點，大量文獻(xiàn)提出了高階FDTD方法來提高計算精度，但高階方法對穩(wěn)定性條件的要求更強(qiáng)；T.Namiki[2]則提出ADI-FDTD來擺脫Courant穩(wěn)定性條件的束縛，但ADI-FDTD方法的數(shù)值色散性較FDTD法差。[3]總而言之，這些方法的效果并不盡如人意，原因在于這些方法破壞了Maxwell方程的辛結(jié)構(gòu)。我們知道Maxwell方程可以看成一個無窮維的Hamilton系統(tǒng)，而Hamilton的算法應(yīng)在辛幾何框架內(nèi)產(chǎn)生，該系統(tǒng)隨時間的演化永遠(yuǎn)是辛變換演進(jìn)，所以正確的離散算法應(yīng)是辛變換，這樣的算法稱為Hamilton算法或辛算法。[4]在對Maxwell方程離散求解時，保持辛結(jié)構(gòu)就顯得非常重要。本文首先詳細(xì)介紹了Hamilton系統(tǒng)中的辛算法，其次將Maxwell方程表述為Hamilton系統(tǒng)，在此基礎(chǔ)上得到保持電磁場辛結(jié)構(gòu)的時域有限差分方法，即辛?xí)r域有限差分（SFDTD）法，最后對其穩(wěn)定性及數(shù)值色散性進(jìn)行了簡要分析。

2 Hamilton系統(tǒng)和辛算法

其中{cl}和{dl}是傳播系數(shù)，m和p分別是辛算法的級數(shù)與階數(shù)。該方法稱之為辛傳播子方法或分解算子方法。為了確定傳播子系數(shù)，將（7）式在t=0點Taylor展開：

比較LHS和RHS兩式中τ的多項式系數(shù)，如果直到τp為止，兩端的系數(shù)都相等，則該辛算法稱之為p階辛算法，由此可確定傳播子系數(shù)。常用的是二級二階、三級三階和五級四階的辛算法，[7,8,9]其傳播子系數(shù)如表1所示。

表1 m級p階辛傳播子系數(shù)

3 Maxwell方程的Hamilton表述

在線性各向同性介質(zhì)中，時域Maxwell方程為：其中，ε，μ分別為相對介電常數(shù)和磁導(dǎo)率。

在Hamilton系統(tǒng)中，方程（10）和（11）可改寫成無窮維Hamilton函數(shù)方程形式，即

4 辛?xí)r域有限差分法（SFDTD）

由以上將Maxwell方程表述為Hamilton函數(shù)，并假設(shè)在時間尺度上，對離散網(wǎng)格點{nΔt，n=1，2，3…}上應(yīng)用（7）式對（13）式近似計算，并根據(jù)（9）式把每個時間間隔 Δt分裂為 m 級，這樣就保持了 Maxwell方程的辛性質(zhì)。在空間方向上的偏微分應(yīng)用離散網(wǎng)格點1/2，…}上的差分進(jìn)行近似計算。

對（13）式，令 LA，LB分別為：

至此，辛格式有效應(yīng)用于Maxwell方程組的時域計算中。

由于（LA）k=0，（LB）k=0（k=2，3…）指數(shù)函數(shù) exp（τLA）和 exp（τLB）可以用下式進(jìn)行數(shù)值計算，

對Maxwell方程組在空間方向上的離散化，將采用差分方法近似空間坐標(biāo)的微分，用jΔy，kΔz；nΔt）表示 f在離散網(wǎng)格點（iΔx， jΔy，kΔz）上第 n 個時間步長上的函數(shù)值。 對空間采用四階差分格式近似偏微分有：

將以上m級p階辛格式和4階空間差分格式相結(jié)合，得到辛?xí)r域有限差分法（SFDTD），[10]記為S（m，p；4），則在x方向上的電磁感應(yīng)強(qiáng)度Dx和Bx的迭代格式可表示為：

CFL的最大值即CFLmax決定了算法穩(wěn)定性的上界，通過優(yōu)化算法，可得到m級p階SFDTD算法的CFLmax列于表2中，作為參考，F(xiàn)DTD方法的CFLmax也同時列出。

另外兩個方向上的電場與磁場分量的迭代格式可類似得到。

5 SFDTD數(shù)值方法的穩(wěn)定性及色散性分析

對SFDTD方法的穩(wěn)定性及數(shù)值色散性的分析，可采用增長矩陣分析方法，[11]為方便分析，在以下的數(shù)值計算過程中，對空間方向上的離散采用均勻網(wǎng)格Δx=Δy=Δz=Δ，定義穩(wěn)定性常數(shù)：

表2 SFDTD算法的CFLmax

用有限差分進(jìn)行差分離散時，將會在離散網(wǎng)格中引起所模擬波模的色散，即在差分網(wǎng)格中，數(shù)值波模的傳播速度將隨頻率而改變，這種色散將導(dǎo)致非物理因素引起的脈沖波形畸變、人為的各向異性及虛假的折射現(xiàn)象。因此，數(shù)值色散是有限差分法中必須考慮的另一個因素。

FDTD算法的數(shù)值色散關(guān)系[1]為：

其中，為角頻率 ω，kx，ky，kz為入射波矢量沿x，y，z方向的分量。

SFDTD算法的數(shù)值色散關(guān)系為：

6 結(jié)論

論文重點闡述了Hamilton系統(tǒng)的辛傳播子算法，給出了m級p階辛算法的傳播子系數(shù)，確保了該離散格式的辛結(jié)構(gòu)的完整性。然后將Maxwell方程組表述為Hamilton格式，為求解該方程組開辟了新的思路，并進(jìn)一步給出了Hamilton系統(tǒng)中Maxwell方程組的時間離散格式，結(jié)合4階空間差分離散格式，得到了電磁場的隨時間和空間推進(jìn)的迭代格式，即辛?xí)r域有限差分方法（SFDTD），最后推導(dǎo)了該方法的穩(wěn)定性條件以及由時間和空間離散引起的數(shù)值色散特性。可以說，SFDTD方法的提出，為今后進(jìn)行電磁散射及輻射問題的分析計算提供了一種精確的數(shù)值方法，可行性很強(qiáng)。

[1]葛德彪，閆玉波.電磁波時域有限差分方法[M].西安：西安電子科技大學(xué)出版社，2001.

[2]T.Namiki.A new FDTD algorithm based on alternating-direction implicit method[J].IEEE Transactions on Microwave Theory and Techniques.1999，47:2003-2007.

[3]F.Zheng，Z.Chen.Numerical dispersion analysis of the unconditionally stable 3D ADI-FDTD method[J].IEEE Transactions on Microwave Theory and Techniques.2001，49:1006-1009.

[4]馮康，秦孟兆.哈密爾頓系統(tǒng)的辛幾何算法[M].杭州：浙江科學(xué)技術(shù)出版社，2003.

[5]Yoshida H.Construction of higher order symplectic integrators[J].Phys Lett A.1990，150:262-268.

[6]M Kusaf，A.Y.Oztoprak，D.S.Daoud.Optimized exponential operator coefficients for symplectic FDTD method[J].IEEE Microwave and Wireless Components Letters，2005，15:86-88.

[7]Kostas T，Simos T E.Symplectic methods for the numerical solution of the radial Shrdinger equation[J].J Chem Phys，2002，34:83-93.

[8]Omelyan I P，Mryglod L M，F(xiàn)olk R.Optimized Forest-Ruth-and Suzuki-like algorithms for interation of motion in manybody systems[J].Comp Phys Comm，2002，146:188-202.

[9]Ruth F D，A canonical integration technique[J].IEEE Trans Nucl Sci，1983，30:2669-2671.

[10]T Hirono，W Lui，S Seki，Y Yoshikuni.A three-dimensional fourth-order finite-difference time-domain scheme using a symplectic integrator propagator[J].IEEE Trans actions on Microwave Theory and Techniques，2001，49:1640-1648.

[11]王秉中.計算電磁學(xué)[M].北京：科學(xué)出版社，2002.

THE NUMERIC SOLUTION OF THE MAXWELL'S EQUATIONS IN THE HAMILTON SYSTEM

LIU Shuang-bing CHEN Hai-bo YANG Han-sheng
（The Department of Physics and Electronics Science，Chaohu College，Chaohu Anhui 238000）

The symplectic algorithm for Hamilton function is expatiated，the Maxwell's equations are written as Hamilton function. Then the idea of keeping conversation of energy and symmetry in the symplectic algorithm is used to solve the Maxwell's equations.Combining with the traditional Finite Difference Time Domain method，we obtain the discrete iterative equations of electromagnetism in time and space domain，named Symplectic Finite Difference Time Domain algorithm.Finally the stability and the numeric dispersion of the SFDTD algorithm are analyzed simply.

Hamilton system；symplectic algorithm；finite difference time domain；numeric dispersion

O441

A

1672－2868（2010）06－0062－05

2010－09-12

巢湖學(xué)院自然科學(xué)研究資助項目（XLY-201010）

劉雙兵（1982-），男，安徽宿松人。碩士研究生，助教，研究方向：電磁場數(shù)值計算。

責(zé)任編輯：宏 彬