本发明涉及一种三维交替迭代无条件稳定fdtd算法，属于fdtd算法技术领域。
背景技术：
时域有限差分算法(finite-differencetime-domain，fdtd)是一种直观的时域电磁场计算方法，目前已经广泛运用到电磁兼容、天线特性、集成电路、吸波材料等领域的数值模拟仿真中。
传统fdtd算法是k.s.yee在1966年提出的，采用电、磁场分量交替网格计算的方法，巧妙地描述了maxwell方程组在空间中的传播。但传统fdtd算法必须满足柯朗(courant-friedrich-lewy，cfl)稳定性条件的要求，当空间步长取值很小时，时间步长也必须相应地取值很小，这大大增加了数值仿真的计算时间，尤其不适合于具有精细网格的计算。
1999年，namiki等提出了基于交替隐式差分方法(alternating-directionimplicitmethod，adi)的fdtd算法。adi-fdtd算法突破了cfl稳定性条件的限制，适用于精细结构的数值仿真，但是数值色散误差的影响较大。2003年，t.k.sarkar等提出了基于laguerre基的fdtd算法，利用laguerre基将电磁场分量展开到laguerre域，并采用galerkin法解耦了时间变量和空间变量的影响。laguerre基fdtd算法不受cfl稳定性条件的限制，是无条件稳定的，数值色散误差也较小，但是该算法需要求解一个大型的线性方程组，需要占用大量的计算机内存，计算效率也受到影响。
技术实现要素：
针对现有技术的不足，本发明的目的在于提供一种三维交替迭代无条件稳定fdtd算法，以解决现有技术中fdtd算法由于受cfl稳定性条件的限制而不适合于具有精细网格的计算、由于需占用大量内存资源而影响计算效率的技术问题。
为解决上述技术问题，本发明所采用的技术方案是：
一种三维交替迭代无条件稳定fdtd算法，其特征是，包括如下步骤：
将三维maxwell方程组变换为一个laguerre域的矩阵形式的maxwell方程；
在所述maxwell方程一侧添加误差项，使所述maxwell方程一侧的系数矩阵分解成两个方阵乘积；
基于系数矩阵分解成两个方阵乘积的maxwell方程求解laguerre域的电场的中间变量，基于所述中间变量求解laguerre域的具有误差的电场分量，基于所述具有误差的电场分量求解laguerre域的具有误差的磁场分量；
利用迭代算法求解所述中间变量、所述具有误差的电场分量、所述具有误差的磁场分量，直至获取所述误差小于预设阈值的laguerre域的电场分量和磁场分量；
基于所述laguerre域的电场分量和磁场分量获取时域的电场分量和磁场分量。
进一步地，将三维maxwell方程组变换为一个laguerre域的矩阵形式的maxwell方程，包括：
利用laguerre正交基将三维maxwell方程组展开到laguerre域，获取laguerre域的六个maxwell方程；
利用构造列向量法将所述laguerre域的六个maxwell方程变换为一个laguerre域的矩阵形式的maxwell方程。
进一步地，基于所述laguerre域的电场分量和磁场分量获取时域的电场分量和磁场分量，包括：利用laguerre正交基展开的逆过程，将所述laguerre域的电场分量和磁场分量变换为时域的电场分量和磁场分量。
进一步地，所述时域的电场分量和磁场分量，其计算公式如下：
式中，ex(r,t)、ey(r,t)、ez(r,t)为时域的电场分量，hx(r,t)、hy(r,t)、hz(r,t)为时域的磁场分量，
进一步地，所述迭代算法为高斯迭代算法。
进一步地，所述laguerre域的矩阵形式的maxwell方程，其表达式如下：
(i-m-n)xq＝yq-1 sq，其中，
式中，(i-m-n)为maxwell方程一侧的系数矩阵，xq为由q阶未知电、磁场分量构成的向量，yq-1为由(q-1)阶以下已知电、磁场分量构成的向量，sq为由q阶场源分量构成的向量，
进一步地，所述误差项，其表达式如下：
mn(xq xq-1)；
式中，xq-1为由(q-1)阶已知电、磁场分量构成的向量。
进一步地，在一侧添加误差项后的laguerre域的矩阵形式的maxwell方程，其表达式如下：
(i-m)(i-n)xq＝-mnxq-1 yq-1 sq；
式中，(i-m)(i-n)为所述两个方阵乘积。
进一步地，所述迭代算法，其计算公式如下：
式中，m为迭代次数，
与现有技术相比，本发明所达到的有益效果：
(1)不受cfl稳定性条件的限制，即无条件稳定。传统的fdtd算法需要满足cfl稳定性条件，时间步长的选取受到空间步长的限制，时间步长选取不当会导致数值计算结果出现重大偏差，甚至无法收敛。特别是对于一些精细结构而言，空间步长须取得很小才能够准确地描述散射体的物理结构，从而导致时间步长也必须取得很小，计算时间大大增加。本发明算法基于laguerre基函数，通过变换将时间变量消除，因而计算时间不受时间步长选取的影响，大大提高了计算效率；
(2)计算效率很高。本发明算法将电、磁场的求解过程解耦，相比于传统的laguerre基fdtd算法，大大降低了矩阵维度。不仅如此，实际求解过程可以发现，每次单独的线性方程组所涉及的未知量全在一条直线上，这样就把三维问题化成一维处理，因而大大减少了计算时间，提高了计算效率；
(3)计算占用内存资源很少。传统的laguerre基fdtd算法由于需要求解大型线性方程组，占用的计算机内存量非常大。本发明算法达到了解耦、降维的效果，因而大大降低了对计算机内存的需求；
(4)具有很高的精度。本发明算法的精度可以通过调整迭代次数进行调节，还可以通过局部迭代，在一定范围内得到更高精度。一般而言，8次整体迭代后就能满足精度要求。
附图说明
图1是本发明实施例中所述窄缝微带线结构示意图；
图2是本发明实施例中采用不同算法得到观测点的电场分量ex的波形；
图3是本发明实施例中不同迭代次数时观测点的电场分量ex的相对误差。
具体实施方式
本发明具体实施方式提供了一种三维交替迭代无条件稳定fdtd算法，包括如下步骤：
步骤1，利用laguerre正交基将三维maxwell方程组展开到laguerre域，利用构造变量的方法将laguerre域的六个maxwell方程写成一个矩阵方程的形式。
在无耗、均匀、各向同性介质中，三维时域maxwell方程组为：
式中，ε为介电常数，μ为磁导率，ex、ey、ez为电场沿x、y、z三个方向的分量，hx、hy、hz为磁场沿x、y、z三个方向的分量，jx、jy、jz为电流源沿x、y、z三个方向的分量。
式(1)～(6)中的电场分量和磁场分量可以用laguerre基函数展开为：
其中，
laguerre基函数具有如下性质：
式中，u(r,t)代表任一电场或磁场分量，δpq为狄拉克函数，up(r)为u(r,t)的p阶laguerre系数，uk(r)为u(r,t)的k阶laguerre系数。因此将式(7)～(12)代入式(1)～(6)，在方程两边同时乘以基函数
其中，
式中，
构造如下变量：
则式(15)～(20)可以写成一个矩阵方程：
(i-m-n)xq＝yq-1 sq(32)
式中，xq为由q阶未知电、磁场分量构成的向量，yq-1为由(q-1)阶以下已知电、磁场分量构成的向量，sq为由q阶场源分量构成的向量，
步骤2，为解耦电场分量和磁场分量的求解过程，在上述矩阵形式的maxwell方程左边主动添加一个引起误差的误差项，使得式(32)左边的系数矩阵可以分解成两个方阵乘积的形式；然后采用分步求解法，将变换后的maxwell方程分解成三步求解，第一步求解电场的中间变量，第二步求解电场，第三部求解磁场。
在式(32)左边加上一个误差项mn(xq xq-1)，并整理得：
(i-m)(i-n)xq＝-mnxq-1 yq-1 sq(33)
式(33)可以分解为两个等式进行计算：
(i-m)x*q＝-nxq-1 yq-1 sq(34)
(i-n)xq＝x*q nxq-1(35)
式中，
利用式(25)～(31)展开式(34)～(35)，可以得到一组求解三维空间电磁场分量的基本表达式：
利用中心差分法可以很方便地求解式(36)～(44)，具体的求解顺序是：首先，求中间变量
需要注意的是，式(36)～(44)中每个式子用中心差分法展开以后都构成一个线性方程组，但方便的是方程左边的未知量只在一个方向上进行差分，所涉及的未知量全在一条直线上，这样就把三维空间的未知量分解成一条直线、一条直线地进行求解，几何级数地减小了系数矩阵维度，提高了计算效率，降低了对计算机内存的需求。
步骤3，步骤2计算出的laguerre域的电场分量
具体做法是：将步骤2的解作为初值
重复以上步骤，可以得到三维fdtd算法的迭代算法：
式中，m为迭代次数。上式可以分解为如下两个等式进行计算：
利用式(25)～(31)，展开式(47)～(48)，得：
式(49)～(57)的求解过程，如同式(36)～(44)。其求解过程实际上是电场中间变量→电场分量→磁场分量→电场中间变量的交替方向迭代方法。为了提高迭代算法的收敛速度，在迭代算法中引入高斯迭代思想，将迭代过程中得到的一些变量提前代入迭代过程，即将式(50)中的
其它方程不变，可以提高迭代效率。当迭代误差达到预设的精度要求时，停止迭代，就可以算出所有电磁场分量的q阶laguerre基函数的系数。最后利用(7)～(12)将信号还原到时域，即可获得计算域内任意场点任意时刻的电磁场波形。
本发明算法的优点在于：
(1)不受cfl稳定性条件的限制，即无条件稳定。传统的fdtd算法需要满足cfl稳定性条件，时间步长的选取受到空间步长的限制，时间步长选取不当会导致数值计算结果出现重大偏差，甚至无法收敛。特别是对于一些精细结构而言，空间步长须取得很小才能够准确地描述散射体的物理结构，从而导致时间步长也必须取得很小，计算时间大大增加。本发明算法基于laguerre基函数，通过变换将时间变量消除，如式(15)～(20)所示，因而计算时间不受时间步长选取的影响，大大提高了计算效率；
(2)计算效率很高。本发明算法将电、磁场的求解过程解耦，相比于传统的laguerre基fdtd算法，大大降低了矩阵维度。不仅如此，实际求解过程可以发现，每次单独的线性方程组所涉及的未知量全在一条直线上，这样就把三维问题化成一维处理，如式(36)～(44)所示，因而大大减少了计算时间，提高了计算效率；
(3)计算占用内存资源很少。传统的laguerre基fdtd算法由于需要求解大型线性方程组，占用的计算机内存量非常大。本发明算法达到了解耦、降维的效果，如式(36)～(44)所示，因而大大降低了对计算机内存的需求；
(4)具有很高的精度。本发明算法的精度可以通过调整迭代次数进行调节，还可以通过局部迭代，在一定范围内得到更高精度。一般而言，8次整体迭代后就能满足精度要求。
下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。
以不连续微带线的数值仿真为例，其结构和尺寸如图1所示，窄缝处采用扩展网格。激励源为高斯脉冲，计算时长为tf＝0.8ns。分别采用传统fdtd算法、adi-fdtd算法、传统laguerre基fdtd算法和本发明算法进行数值仿真。
表1给出了不同计算方案重要参数设置和计算结果。其中ns、nt分别表示整体、局部迭代次数，局部迭代的区域为导带断点前后的4×10×46个网格。从表1中可以看到，本发明算法在ns＝3、nt＝9时，计算时间仅为6.97s，是传统fdtd算法的8.9％，是laguerre基fdtd算法的23.8％。从内存需求上看，本发明算法为2.82mb，且不随迭代次数而改变，比传统fdtd算法略大，是laguerre基fdtd算法的8.6％。
表1：不同算法消耗的计算资源
图2为采用不同方法计算的观测点的电场分量ex的波形，可以看出，本发明算法在ns＝3、nt＝5时，波形就与其它三种算法吻合的很好。
本专利算法的误差可以通过相对误差来判定：
式中，rdb为以分贝为单位定义的相对误差，
本发明算法结合了adi算法和利用laguerre基展开的思想，将laguerre域的maxwell方程组写成矩阵形式，通过主动增加误差项的方法，解耦了电、磁场分量的计算。更加有利的是，这种方法得到线性方程组中的未知场量全在一条直线上，即将需要求解的三维空间场问题化简到一条直线、一条直线上进行求解，几何级数地降低了线性方程组的维度，大大提高了计算效率，减小了内存消耗。最后通过交替方向迭代的方法，使误差项的大小趋于零，得到电磁场分量的精确解。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。
1.一种三维交替迭代无条件稳定fdtd算法，其特征是，包括如下步骤：
将三维maxwell方程组变换为一个laguerre域的矩阵形式的maxwell方程；
在所述maxwell方程一侧添加误差项，使所述maxwell方程一侧的系数矩阵分解成两个方阵乘积；
基于系数矩阵分解成两个方阵乘积的maxwell方程求解laguerre域的电场的中间变量，基于所述中间变量求解laguerre域的具有误差的电场分量，基于所述具有误差的电场分量求解laguerre域的具有误差的磁场分量；
利用迭代算法求解所述中间变量、所述具有误差的电场分量、所述具有误差的磁场分量，直至获取所述误差小于预设阈值的laguerre域的电场分量和磁场分量；
基于所述laguerre域的电场分量和磁场分量获取时域的电场分量和磁场分量。
2.根据权利要求1所述的三维交替迭代无条件稳定fdtd算法，其特征是，将三维maxwell方程组变换为一个laguerre域的矩阵形式的maxwell方程，包括：
利用laguerre正交基将三维maxwell方程组展开到laguerre域，获取laguerre域的六个maxwell方程；
利用构造列向量法将所述laguerre域的六个maxwell方程变换为一个laguerre域的矩阵形式的maxwell方程。
3.根据权利要求2所述的三维交替迭代无条件稳定fdtd算法，其特征是，基于所述laguerre域的电场分量和磁场分量获取时域的电场分量和磁场分量，包括：利用laguerre正交基展开的逆过程，将所述laguerre域的电场分量和磁场分量变换为时域的电场分量和磁场分量。
4.根据权利要求3所述的三维交替迭代无条件稳定fdtd算法，其特征是，所述时域的电场分量和磁场分量，其计算公式如下：
式中，ex(r,t)、ey(r,t)、ez(r,t)为时域的电场分量，hx(r,t)、hy(r,t)、hz(r,t)为时域的磁场分量，
5.根据权利要求1所述的三维交替迭代无条件稳定fdtd算法，其特征是，所述迭代算法为高斯迭代算法。
6.根据权利要求1所述的三维交替迭代无条件稳定fdtd算法，其特征是，所述laguerre域的矩阵形式的maxwell方程，其表达式如下：
(i-m-n)xq＝yq-1 sq，其中，
式中，(i-m-n)为maxwell方程一侧的系数矩阵，xq为由q阶未知电、磁场分量构成的向量，yq-1为由(q-1)阶以下已知电、磁场分量构成的向量，sq为由q阶场源分量构成的向量，
7.根据权利要求6所述的三维交替迭代无条件稳定fdtd算法，其特征是，所述误差项，其表达式如下：
mn(xq xq-1)；
式中，xq-1为由(q-1)阶已知电、磁场分量构成的向量。
8.根据权利要求7所述的三维交替迭代无条件稳定fdtd算法，其特征是，在一侧添加误差项后的laguerre域的矩阵形式的maxwell方程，其表达式如下：
(i-m)(i-n)xq＝-mnxq-1 yq-1 sq；
式中，(i-m)(i-n)为所述两个方阵乘积。
9.根据权利要求8所述的三维交替迭代无条件稳定fdtd算法，其特征是，所述迭代算法，其计算公式如下：
式中，m为迭代次数，
