CN105825015A - 一种用于磁化等离子体的时域有限差分方法 - Google Patents

Abstract

本发明公开了一种用于磁化等离子体的时域有限差分方法，具体为：步骤1、输入模型文件；步骤2、初始化和设置步骤1中的参数；步骤3、利用步骤2的参数，计算电场分量系数步骤4、利用步骤2的参数，计算电场分量系数步骤5、添加场源到y方向上的磁场分量系数中，并利用步骤3所得计算磁场分量系数步骤6、利用步骤4所得计算步骤7、利用电场分量系数计算极化电流密度步骤8、更新计算电磁场分量系数的辅助变量；步骤9、更新计算观测点处的电磁场分量；步骤10、将q+1赋值给q，并判断q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。本发明计算速度快，内存消耗小，且对于低频和凋落波具有很好的吸收效果。

Description

一种用于磁化等离子体的时域有限差分方法

技术领域

本发明属于计算电磁学技术领域，具体涉及一种用于磁化等离子体的时域有限差分方法。

背景技术

时域有限差分(Finite-differencetime-domain，FDTD)方法广泛用于求解磁化等离子体中的电磁问题，但是由于其时间步长受柯西稳定性条件的限制，不能取的太大，方法时效性不好，计算速度较慢。为了消除柯西稳定性条件的限制，人们提出了无条件稳定时域有限差分方法，比如：交替方向隐式(Alternating-Direction-Implicit，ADI)的时域有限差分(ADI-FDTD)方法和基于加权拉盖尔多项式的时域有限差分(Finite-differencetime-domainwithWeighted-Laguerre-polynomials，WLP-FDTD)方法。在这些方法中，ADI-FDTD方法在使用较大的时间步长时会产生很大的色散误差，而WLP-FDTD方法既能消除柯西稳定性条件的限制，又能解决ADI-FDTD方法在使用较大的时间步长时会产生很大的色散误差这个难题，因此WLP-FDTD方法可以高效的求解磁化等离子体中的电磁问题，然而，这种WLP-FDTD方法在求解电磁场过程中，会产生一个大型的稀疏矩阵方程，直接求解此方程会使得计算较复杂，内存消耗较大，于是提出了一种因式分解的WLP-FDTD方法。

另外由于计算机容量的限制，电磁场的计算只能在有限区域进行。为了能模拟开域电磁波传播过程，必须在计算区域的截断边界处给出吸收边界条件。有人提出了完全匹配层(Perfectlymatchedlayer，PML)吸收边界，后来PML被广泛应用于计算区域的截断，而且被证明是非常有效的，但是研究发现这种传统PML对低频以及凋落波的吸收效果并不理想；使用带有复频率偏移(Complexfrequencyshift，CFS)因子的扩展坐标的PML吸收边界可以有效地改善传统PML对低频，凋落波的吸收效果。

发明内容

本发明的目的是提供一种用于磁化等离子体的时域有限差分方法，解决了现有求解磁化等离子体中的电磁问题时存在的计算速度慢、对低频以及凋落波吸收效果不好的问题。

本发明所采用的技术方案是，一种用于磁化等离子体的时域有限差分方法，具体按照以下步骤：

步骤1、输入模型文件；

步骤2、初始化和设置步骤1内模型文件的参数；

步骤3、利用步骤2的参数，更新计算整个计算区域x方向上电场分量系数

步骤4、利用步骤2的参数，更新计算整个计算区域y方向上电场分量系数

步骤5、添加场源到y方向上的磁场分量系数中，并利用步骤3所得x方向上电场分量系数更新计算整个计算区域y方向上的磁场分量系数

步骤6、利用步骤4所得y方向上电场分量系数更新计算整个计算区域的x方向上磁场分量系数

步骤7、利用步骤3和步骤4所得电场分量系数，更新计算整个计算区域的极化电流密度

步骤8、更新计算整个计算区域的电磁场分量系数的辅助变量，其中电磁场分量包括E_x,E_y,H_x,H_y；

步骤9、更新计算观测点处的电磁场分量，电磁场分量包括：E_x,E_y,H_x,H_y；

步骤10、将q+1赋值给步骤9所得电磁场分量中的q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。

本发明的特点还在于：

步骤1输入模型文件，具体输入的参数包括：

计算区域大小N_z；空间步长Δz；时间步长Δt；真空中的电导率σ、磁导率μ₀、介电常数ε₀；等离子体中的碰撞频率ν；等离子体频率ω_p；电子回旋频率ω_b；等离子体在计算区域中的位置；吸收边界层数NPML与相关参数κ_zmax，α_zmax，σ_zmax，其中，κ_zmax取整数，κ_zmax取值范围为[1，60]，α_zmax取值范围为[0，1)，σ_zmax/σ_opt取值范围为(0，12]，σ_opt＝(m+1)/150πΔz，m取值范围为[1，20]，Δz取值范围为λ为源的波长；仿真计算时长T_f；加权拉盖尔多项式的阶数q(q≥0且为整数)；时间尺度因子s，其中s取值范围为[10⁹,10¹³]；观测点；场源参数。

步骤2初始化参数，具体为：

将整个计算区域的电磁场分量系数整个计算区域的极化电流密度整个计算区域的电磁场分量系数的和整个计算区域的极化电流密度的和整个计算区域的辅助变量其中F_ζ表示E_x,E_y,H_x,H_y、和拉盖尔多项式全部初始化为零、等离子体参数(p₁,p₂,p₃)初始化为p₁＝0，p₂＝0，PML系数(C_1z,C_2z,C₃)初始化为C_1z＝1/(1+0.5ε₀s)，C_2z＝1，C₃＝ε₀/μ₀；

设置参数具体为：

设置带有CFS因子的SC-PML吸收边界的参数σ_z,κ_z,α_z：

σ_z＝σ_zmax|z-z₀|^m/d^m

κ_z＝1+(κ_zmax-1)|z-z₀|^m/d^m

α_z＝α_zmax

式中z₀为PML层与非PML截面位置，d是PML吸收边界的厚度；

设置PML系数C_1z,C_2z：

C_1z＝1/(κ_zα_z+σ_z+0.5κ_zε₀s)

C_2z＝(2α_z/ε₀s+1)

设置等离子体参数：

$p_1=\frac{{\mathrm{\ω}}_b}{0.5s+\mathrm{\υ}},p_2=\frac{4{\mathrm{\ω}}_p^2}{(s^2+2s\mathrm{\υ})(1+p_1^2)},p_3=\frac{4}{{\mathrm{\ϵ}}_0(s+2\mathrm{\υ})(1+p_1^2)}.$

步骤3利用步骤2的参数，更新计算整个计算区域x方向上电场分量系数具体为：

步骤3.1、将步骤2的参数代入电场分量系数在计算区域的方程，具体为：

$\begin{array}{c}-\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}{E}_y^q{|}_{k-1}\\ +(1+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k})E_y^q{|}_k\\ -\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k}{E}_y^q{|}_{k+1}\\ =\left(\begin{array}{c}{V}_{Ey}^{q-1}{|}_k-{\mathrm{bV}}_{Ex}^{q-1}{|}_k+\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}({\mathrm{bV}}_{Hy}^{q-1}{|}_{k+1/2}-{\mathrm{bV}}_{Hy}^{q-1}{|}_{k-1/2})\\ +\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}(V_{Hx}^{q-1}{|}_{k+1/2}-V_{Hx}^{q-1}{|}_{k-1/2})\end{array}\right)/(1+b^2)\end{array}$

式中，k表示第k个计算网格，

步骤3.2、使用追赶法求解步骤3.1的方程，得到整个计算区域x方向上的电场分量系数

步骤4利用步骤2的参数，更新计算整个计算区域y方向上电场分量系数具体为：

步骤4.1、将步骤2的参数代入电场分量系数在计算区域的方程，具体为：

$\begin{array}{c}-\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}{E}_x^q{|}_{k-1}\\ +(1+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k})E_x^q{|}_k\\ -\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k}{E}_x^q{|}_{k+1}\\ =\left(\begin{array}{c}{V}_{Ex}^{q-1}{|}_k+{\mathrm{bV}}_{Ey}^{q-1}{|}_k+\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}({\mathrm{bV}}_{Hx}^{q-1}{|}_{k+1/2}-{\mathrm{bV}}_{Hx}^{q-1}{|}_{k-1/2})\\ +\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}(V_{Hy}^{q-1}{|}_{k+1/2}-V_{Hy}^{q-1}{|}_{k-1/2})\end{array}\right)/(1+b^2)\end{array}$

步骤4.2、使用追赶法求解步骤4.1的方程，得到整个计算区域y方向上的电场分量系数

步骤5添加场源到y方向上的磁场分量系数中，并利用步骤3所得x方向上电场分量系数更新计算整个计算区域的y方向上磁场分量系数具体为：

场源的表达式为：

I_my(t)＝exp(-(t-t₀)²/τ²)

式中，t₀,τ为场源参数；

将步骤3所得x方向上电场分量系数代入磁场分量系数计算公式中，更新磁场分量系数

$H_y^q{|}_{k+1/2}=V_{Hy}^{q-1}{|}_{k+1/2}-C_3\frac{C_{2z}{|}_{k+1/2}}{{\mathrm{\Δz}}_{k+1/2}}(E_x^q{|}_{k+1}-E_x^q{|}_k)$

式中

步骤6利用步骤4所得y方向上电场分量系数更新计算整个计算区域的x方向上磁场分量系数具体为

$H_x^q{|}_{k+1/2}=V_{Hx}^{q-1}{|}_{k+1/2}+C_3\frac{C_{2z}{|}_{k+1/2}}{{\mathrm{\Δz}}_{k+1/2}}(E_y^q{|}_{k+1}-E_y^q{|}_k)$

式中

步骤7利用步骤3和步骤4所得电场分量系数，更新计算整个计算区域的极化电流密度具体更新公式为：

$J_x^q{|}_k==0.5{\mathrm{\ϵ}}_{0}s(p_2{E}_x^q{|}_k-p_1{p}_2{E}_y^q{|}_k-p_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_x^k{|}_k+p_1{p}_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_y^k{|}_k)$

$J_y^q{|}_k=0.5{\mathrm{\ϵ}}_{0}s(p_2{E}_y^q{|}_k+p_1{p}_2{E}_x^q{|}_k-p_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_y^k{|}_k-p_1{p}_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_x^k{|}_k).$

步骤8更新计算整个计算区域的电磁场分量系数的辅助变量，具体更新公式为：

步骤9更新计算观测点处的电磁场分量，具体按照以下公式更新计算：

式中U表示电磁场分量E_x,E_y,H_x,H_y，U^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s＞0的扩展时间，是q阶拉盖尔多项式。

本发明的有益效果是：

1)在直角坐标系下，通过用加权拉盖尔多项式表示电磁场分量，来解时域麦克斯韦方程，使得在更新计算整个计算区域的电磁场分量系数时不涉及到时间步长，只是在最后计算观测点处的电磁场分量时用到时间步长，因此计算过程中时间步长可以取得比柯西稳定性条件限制的时间步长更大；

2)在求解电磁场分量系数时，将大型稀疏矩阵方程分裂成两个三对角矩阵方程，使得它在计算时比WLP-FDTD方法更简单、计算速度更快、内存消耗更少而且可以对大区域的电磁场问题进行求解；

3)在设置PML系数时，由于采用了CFS因子，并且通过调整CFS因子中的参数，可以使得该吸收边界对低频与凋落波的吸收更加有效；

4)由于采用了复扩展坐标系，使得PML在实现时避免了场的分裂且与媒质无关。

附图说明

图1是本发明的流程示意图；

图2是本发明实施例中本发明方法、解析解、WLP-FDTD方法计算的左旋极化波的反射系数振幅图；

图3是本发明实施例中本发明方法与解析解、WLP-FDTD方法计算的左旋极化波的透射系数振幅图；

图4是本发明方法与WLP-FDTD方法CPU时间随网格大小变化的图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

本发明一种用于磁化等离子体的时域有限差分方法，磁化等离子体中，电磁场所满足的复扩展坐标系下的麦克斯韦方程推导过程如下：

在各向异性色散介质碰撞磁化等离子体中，麦克斯韦方程组和相关的联立方程为

$\▿\×H={\mathrm{\ϵ}}_0\frac{\∂E}{\∂t}+J---\left(1\right)$

$\▿\×E=-{\mathrm{\μ}}_0\frac{\∂H}{\∂t}---\left(2\right)$

$\frac{\∂J}{\∂t}+\mathrm{\υ}J={\mathrm{\ϵ}}_0{\mathrm{\ω}}_p^{2}E+{\mathrm{\ω}}_b\×J---\left(3\right)$

式中，H是磁场强度，E是电场强度，J是极化电流密度，ε₀、μ₀分别为真空中的介电常数和磁导率，υ是等离子体碰撞频率，是等离子体角频率的平方，ω_b＝eB₀/m是电子回旋频率，B₀是外磁场，e、m分别是电子的电量和质量。

设外磁场的方向为+z轴，方程(3)可写为

$\frac{\∂J_x}{\∂t}+{\mathrm{\υJ}}_x={\mathrm{\ϵ}}_0{\mathrm{\ω}}_p^2{E}_x-{\mathrm{\ω}}_b{J}_y---\left(4\right)$

$\frac{\∂J_y}{\∂t}+{\mathrm{\υJ}}_y={\mathrm{\ϵ}}_0{\mathrm{\ω}}_p^2{E}_y+{\mathrm{\ω}}_b{J}_x---\left(5\right)$

应用扩展坐标的CFS-PML，仅考虑二维TEz的情况，上述麦克斯韦方程组和相关的联立方程可化为：

$-\frac{1}{s_z}\frac{\∂H_y}{\∂z}={\mathrm{\ϵ}}_0\frac{\∂E_x}{\∂t}+J_x---\left(6\right)$

$\frac{1}{s_z}\frac{\∂H_x}{\∂z}={\mathrm{\ϵ}}_0\frac{\∂E_y}{\∂t}+J_y---\left(7\right)$

${\mathrm{\μ}}_0\frac{\∂H_x}{\∂t}=\frac{1}{s_z}\frac{\∂E_y}{\∂z}---\left(8\right)$

${\mathrm{\μ}}_0\frac{\∂H_y}{\∂t}=-\frac{1}{s_z}\frac{\∂E_x}{\∂z}---\left(9\right)$

$\frac{\∂J_x}{\∂t}+{\mathrm{\υJ}}_x={\mathrm{\ϵ}}_0{\mathrm{\ω}}_p^2{E}_x-{\mathrm{\ω}}_b{J}_y---\left(10\right)$

$\frac{\∂J_y}{\∂t}+{\mathrm{\υJ}}_y={\mathrm{\ϵ}}_0{\mathrm{\ω}}_p^2{E}_y+{\mathrm{\ω}}_b{J}_x---\left(11\right)$

式中s_z是带有CFS因子的坐标伸缩因子，取为：

下面我们引入几个辅助变量如下：

${\tilde{H}}_{yz}=\frac{1}{s_z}\frac{\∂H_y}{\∂z}---\left(13\right)$

${\tilde{H}}_{xz}=\frac{1}{s_z}\frac{\∂H_x}{\∂z}---\left(14\right)$

${\tilde{E}}_{xz}=\frac{1}{s_z}\frac{\∂E_x}{\∂z}---\left(15\right)$

${\tilde{E}}_{yz}=\frac{1}{s_z}\frac{\∂E_y}{\∂z}---\left(16\right)$

将式(12)分别代入(13)-(16)，然后做如下变化可得：

$({\mathrm{\κ}}_z{\mathrm{\α}}_z+{\mathrm{\σ}}_z){\tilde{H}}_{yz}+{\mathrm{\κ}}_z{\mathrm{\ϵ}}_0\frac{\∂{\tilde{H}}_{yz}}{\∂t}={\mathrm{\α}}_z\frac{\∂H_y}{\∂z}+{\mathrm{\ϵ}}_0\frac{\∂}{\∂t}\left(\frac{\∂H_y}{\∂z}\right)---\left(17\right)$

$({\mathrm{\κ}}_z{\mathrm{\α}}_z+{\mathrm{\σ}}_z){\tilde{H}}_{xz}+{\mathrm{\κ}}_z{\mathrm{\ϵ}}_0\frac{\∂{\tilde{H}}_{xz}}{\∂t}={\mathrm{\α}}_z\frac{\∂H_x}{\∂z}+{\mathrm{\ϵ}}_0\frac{\∂}{\∂t}\left(\frac{\∂H_x}{\∂z}\right)---\left(18\right)$

$({\mathrm{\κ}}_z{\mathrm{\α}}_z+{\mathrm{\σ}}_z){\tilde{E}}_{xz}+{\mathrm{\κ}}_z{\mathrm{\ϵ}}_0\frac{\∂{\tilde{E}}_{xz}}{\∂t}={\mathrm{\α}}_z\frac{\∂E_x}{\∂z}+{\mathrm{\ϵ}}_0\frac{\∂}{\∂t}\left(\frac{\∂E_x}{\∂z}\right)---\left(19\right)$

$({\mathrm{\κ}}_z{\mathrm{\α}}_z+{\mathrm{\σ}}_z){\tilde{E}}_{yz}+{\mathrm{\κ}}_z{\mathrm{\ϵ}}_0\frac{\∂{\tilde{E}}_{yz}}{\∂t}={\mathrm{\α}}_z\frac{\∂E_y}{\∂z}+{\mathrm{\ϵ}}_0\frac{\∂}{\∂t}\left(\frac{\∂E_y}{\∂z}\right)---\left(20\right)$

我们可以知道

式中U代表E_x、E_y、H_x、H_y，是加权拉盖尔多项式， q≥0；t≥0是q阶拉盖尔多项式。

将式(21)、式(22)代入(17)-(20)，应用加勒金测试过程，得：

式中

$C_{1z}=\frac{1}{({\mathrm{\κ}}_z{\mathrm{\α}}_z+{\mathrm{\σ}}_z+0.5{\mathrm{\κ}}_z{\mathrm{\ϵ}}_{0}s)}---\left(31\right)$

式中，s＞0是时间尺度因子，q是加权拉盖尔多项式的阶数。

将式(21)、式(22)代入(6)-(11)，再应用加勒金的测试过程得到：

$J_x^q=\frac{{\mathrm{\ϵ}}_0{\mathrm{\ω}}_p^2{E}_x^q}{0.5s+\mathrm{\υ}}-\frac{{\mathrm{\ω}}_b}{0.5s+\mathrm{\υ}}{J}_y^q-\frac{s}{0.5s+\mathrm{\υ}}\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_x^k---\left(36\right)$

$J_y^q=\frac{{\mathrm{\ϵ}}_0{\mathrm{\ω}}_p^2{E}_y^q}{0.5s+\mathrm{\υ}}+\frac{{\mathrm{\ω}}_b}{0.5s+\mathrm{\υ}}{J}_x^q-\frac{s}{0.5s+\mathrm{\υ}}\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_y^k---\left(37\right)$

式中，

将(36)代入(37)得到

$J_y^q=0.5{\mathrm{\ϵ}}_{0}s(p_2{E}_y^q+p_1{p}_2{E}_x^q-p_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_y^k-p_1{p}_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_x^k)---\left(38\right)$

将(37)代入(36)得到

$J_x^q==0.5{\mathrm{\ϵ}}_{0}s(p_2{E}_x^q-p_1{p}_2{E}_y^q-p_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_x^k+p_1{p}_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_y^k)---\left(39\right)$

将(35)和(39)代入(32)得到

将(38)和(34)代入(33)得到

将方程(40)、(41)、(34)和(35)放置一起

整理后得到

$E_x^q-\frac{p_1{p}_2}{1+p_2}{E}_y^q+\frac{C_{2z}{D}_z}{1+p_2}{H}_y^q=V_{Ex}^{q-1}$

$E_y^q+\frac{p_1{p}_2}{1+p_2}{E}_x^q-\frac{C_{2z}{D}_z}{1+p_2}{H}_x^q=V_{Ey}^{q-1}---\left(43\right)$

$H_x^q-C_3{C}_{2z}{D}_z{E}_y^q=V_{Hx}^{q-1}$

$H_y^q+C_3{C}_{2z}{D}_z{E}_x^q=V_{Hy}^{q-1}$

其中

将式(43)写成矩阵的形式

(I-A-B)X^q＝V^q-1(45)

式中

$B=\left[\begin{array}{cccc}& b& & \\ -b& & & \\ & & & \\ & & & \end{array}\right],A=\left[\begin{array}{cccc}& & & -a_1\\ & & a_1& \\ & a_2& & \\ -a_2& & & \end{array}\right]---\left(46\right)$

$X^q={\left[\begin{array}{cccc}{E}_x^q& E_y^q& H_x^q& H_y^q\end{array}\right]}^T,V^{q-1}={\left[\begin{array}{cccc}{V}_{Ex}^{q-1}& V_{Ey}^{q-1}& V_{Hx}^{q-1}& V_{Hy}^{q-1}\end{array}\right]}^T$

b＝p₁p₂/(1+p₂),a₁＝C_2zD_z/(1+p₂),a₂＝C_2zC₃D_z

添加微扰项BA(X^q-V^q-1)到(45)得到

(I-B)(I-A)X^q＝V^q-1+BAV^q-1(47)

引进中间变量Y^q＝(I-A)X^q+AV^q-1，于是容易得到

(I-B)Y^q＝(I+A)V^q-1(48)

将上式中的前两个方程展开得到

$\begin{array}{c}{\stackrel{\‾}{E}}_x^q-b{\stackrel{\‾}{E}}_y^q=V_{Ex}^{q-1}-a_1{V}_{Hy}^{q-1}\\ b{\stackrel{\‾}{E}}_x^q+{\stackrel{\‾}{E}}_y^q=V_{Ey}^{q-1}+a_1{V}_{Hx}^{q-1}\end{array}---\left(49\right)$

${\stackrel{\‾}{H}}_x^q=a_2{V}_{Ey}^{q-1}+V_{Hx}^{q-1}$

${\stackrel{\‾}{H}}_y^q=-a_2{V}_{Ex}^{q-1}+V_{Hy}^{q-1}---\left(50\right)$

上式解得

${\stackrel{\‾}{E}}_x^q=b{\stackrel{\‾}{E}}_y^q+V_{Ex}^{q-1}-a_1{V}_{Hy}^{q-1}$

${\stackrel{\‾}{E}}_y^q=-b{\stackrel{\‾}{E}}_x^q+V_{Ey}^{q-1}+a_1{V}_{Hx}^{q-1}$

${\stackrel{\‾}{E}}_x^q=b(-b{\stackrel{\‾}{E}}_x^q+V_{Ey}^{q-1}+a_1{V}_{Hx}^{q-1})+V_{Ex}^{q-1}-a_1{V}_{Hy}^{q-1}$

${\stackrel{\‾}{E}}_y^q=-b(b{\stackrel{\‾}{E}}_y^q+V_{Ex}^{q-1}-a_1{V}_{Hy}^{q-1})+V_{Ey}^{q-1}+a_1{V}_{Hx}^{q-1}---\left(51\right)$

${\stackrel{\‾}{E}}_x^q=(V_{Ex}^{q-1}+{\mathrm{bV}}_{Ey}^{q-1}+a_1{\mathrm{bV}}_{Hx}^{q-1}-a_1{V}_{Hy}^{q-1})/(1+b^2)$

${\stackrel{\‾}{E}}_y^q=(V_{Ey}^{q-1}-{\mathrm{bV}}_{Ex}^{q-1}+a_1{\mathrm{bV}}_{Hy}^{q-1}+a_1{V}_{Hx}^{q-1})/(1+b^2)$

${\stackrel{\‾}{H}}_x^q=a_2{V}_{Ey}^{q-1}+V_{Hx}^{q-1}$

${\stackrel{\‾}{H}}_y^q=-a_2{V}_{Ex}^{q-1}+V_{Hy}^{q-1}---\left(52\right)$

将Y^q＝(I-A)X^q+AV^q-1展开得到

${\stackrel{\‾}{E}}_x^q=E_x^q+a_1{H}_y^q-a_1{V}_{Hy}^{q-1}$

${\stackrel{\‾}{E}}_y^q=E_y^q-a_1{H}_x^q+a_1{V}_{Hx}^{q-1}$

${\stackrel{\‾}{H}}_x^q=-a_2{E}_y^q+H_x^q+a_2{V}_{Ey}^{q-1}$

${\stackrel{\‾}{H}}_y^q=a_2{E}_x^q+H_y^q-a_2{V}_{Ex}^{q-1}---\left(53\right)$

将(51)、(52)和(53)联合解得

$E_x^q-a_1{a}_2{E}_x^q=(V_{Ex}^{q-1}+{\mathrm{bV}}_{Ey}^{q-1}+a_1{\mathrm{bV}}_{Hx}^{q-1}-a_1{V}_{Hy}^{q-1})/(1+b^2)$

$E_y^q-a_1{a}_2{E}_y^q=(V_{Ey}^{q-1}-{\mathrm{bV}}_{Ex}^{q-1}+a_1{\mathrm{bV}}_{Hy}^{q-1}+a_1{V}_{Hx}^{q-1})/(1+b^2)---\left(54\right)$

$H_x^q=V_{Hx}^{q-1}+a_2{E}_y^q$

$H_y^q=V_{Hy}^{q-1}-a_2{E}_x^q$

将b＝p₁p₂/(1+p₂),a₁＝C_2zD_z/(1+p₂),a₂＝C_2zC₃D_z代入上式然后进行中心差分得到

$\begin{array}{c}-\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}{E}_x^q{|}_{k-1}\\ +(1+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k})E_x^q{|}_k\\ -\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k}{E}_x^q{|}_{k+1}\\ =\left(\begin{array}{c}{V}_{Ex}^{q-1}{|}_k+{\mathrm{bV}}_{Ey}^{q-1}{|}_k+\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}({\mathrm{bV}}_{Hx}^{q-1}{|}_{k+1/2}-{\mathrm{bV}}_{Hx}^{q-1}{|}_{k-1/2})\\ +\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}(V_{Hy}^{q-1}{|}_{k+1/2}-V_{Hy}^{q-1}{|}_{k-1/2})\end{array}\right)/(1+b^2)\end{array}---\left(55\right)$

$\begin{array}{c}-\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}{E}_y^q{|}_{k-1}\\ +(1+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k})E_y^q{|}_k\\ -\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k}{E}_y^q{|}_{k+1}\\ =\left(\begin{array}{c}{V}_{Ey}^{q-1}{|}_k-{\mathrm{bV}}_{Ex}^{q-1}{|}_k+\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}({\mathrm{bV}}_{Hy}^{q-1}{|}_{k+1/2}-{\mathrm{bV}}_{Hy}^{q-1}{|}_{k-1/2})\\ +\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}(V_{Hx}^{q-1}{|}_{k+1/2}-V_{Hx}^{q-1}{|}_{k-1/2})\end{array}\right)/(1+b^2)\end{array}---\left(56\right)$

$\begin{array}{c}{H}_x^q{|}_{k+1/2}=V_{Hx}^{q-1}{|}_{k+1/2}+C_3\frac{C_{2z}{|}_{k+1/2}}{{\mathrm{\Δz}}_{k+1/2}}(E_y^q{|}_{k+1}-E_y^q{|}_k)\\ H_y^q{|}_{k+1/2}=V_{Hy}^{q-1}{|}_{k+1/2}-C_3\frac{C_{2z}{|}_{k+1/2}}{{\mathrm{\Δz}}_{k+1/2}}(E_x^q{|}_{k+1}-E_x^q{|}_k)\end{array}---\left(57\right)$

上面三式中，k表示第k个计算网格；在整个计算区域上，式(55)和式(56)可以写成三对角矩阵差分方程，与WLP-FDTD方法相比，这种因式分解的WLP-FDTD方法将大型稀疏矩阵方程的求解转变成两个三对角矩阵方程的求解，于是可以使用追赶法，非常简单的解得整个计算区域电磁场分量系数，最后通过式(21)解得观测点的电磁场分量。

本发明一种用于磁化等离子体的时域有限差分方法，流程如图1所示，具体按照以下步骤：

步骤1输入模型文件，具体输入的参数包括：

计算区域大小N_z；空间步长Δz；时间步长Δt；真空中的电导率σ、磁导率μ₀、介电常数ε₀；等离子体中的碰撞频率υ；等离子体频率ω_p；电子回旋频率ω_b；等离子体在计算区域中的位置；吸收边界层数NPML与相关参数κ_zmax，α_zmax，σ_zmax，其中，κ_zmax取整数，κ_zmax取值范围为[1，60]，α_zmax取值范围为[0，1)，σ_zmax/σ_opt取值范围为(0，12]，σ_opt＝(m+1)/150πΔz，m取值范围为[1，20]，Δz取值范围为λ为源的波长；仿真计算时长Tf；加权拉盖尔多项式的阶数q(q≥0且为整数)；时间尺度因子s，其中s取值范围为[10⁹,10¹³]；观测点；场源参数。

步骤2初始化和设置参数，具体为：

将整个计算区域的电磁场分量系数整个计算区域的极化电流密度整个计算区域的电磁场分量系数的和整个计算区域的极化电流密度的和整个计算区域的辅助变量其中F_ζ表示E_x,E_y,H_x,H_y、和拉盖尔多项式全部初始化为零、等离子体参数(p₁,p₂,p₃)初始化为p₁＝0，p₂＝0，PML系数(C_1z,C_2z,C₃)初始化为C_1z＝1/(1+0.5ε₀s)，C_2z＝1，C₃＝ε₀/μ₀；

设置参数具体为：

设置带有CFS因子的SC-PML吸收边界的参数σ_z,κ_z,α_z：

σ_z＝σ_zmax|z-z₀|^m/d^m

κ_z＝1+(κ_zmax-1)|z-z₀|^m/d^m

α_z＝α_zmax

式中z₀为PML层与非PML截面位置，d是PML吸收边界的厚度，；

设置PML系数C_1z,C_2z：

C_1z＝1/(κ_zα_z+σ_z+0.5κ_zε₀s)

C_2z＝(2α_z/ε₀s+1)

设置等离子体参数：

$p_1=\frac{{\mathrm{\ω}}_b}{0.5s+\mathrm{\υ}},p_2=\frac{4{\mathrm{\ω}}_p^2}{(s^2+2s\mathrm{\υ})(1+p_1^2)},p_3=\frac{4}{{\mathrm{\ϵ}}_0(s+2\mathrm{\υ})(1+p_1^2)}.$

步骤3利用步骤2的参数，更新计算整个计算区域的x方向上电场分量系数，具体为：

步骤3.1、将步骤2的参数代入电场分量系数在计算区域的方程，具体为：

$\begin{array}{c}-\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}{E}_y^q{|}_{k-1}\\ +(1+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k})E_y^q{|}_k\\ -\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k}{E}_y^q{|}_{k+1}\\ =\left(\begin{array}{c}{V}_{Ey}^{q-1}{|}_k+{\mathrm{bV}}_{Ex}^{q-1}{|}_k+\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}({\mathrm{bV}}_{Hy}^{q-1}{|}_{k+1/2}-{\mathrm{bV}}_{Hy}^{q-1}{|}_{k-1/2})\\ +\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}(V_{Hy}^{q-1}{|}_{k+1/2}-V_{Hx}^{q-1}{|}_{k-1/2})\end{array}\right)/(1+b^2)\end{array}$

式中，k表示第k个计算网格，

步骤3.2、使用追赶法求解步骤3.1的方程，得到整个计算区域x方向上的电场分量系数

步骤4利用步骤2的参数，更新计算整个计算区域的y方向上电场分量系数具体为：

步骤4.1、将步骤2的参数代入电场分量系数在计算区域的方程，具体为：

$\begin{array}{c}-\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}{E}_x^q{|}_{k-1}\\ +(1+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k-1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k-1/2}{\mathrm{\Δz}}_k}+\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k})E_x^q{|}_k\\ -\frac{C_3{C}_{2z}{|}_k{C}_{2z}{|}_{k+1/2}}{(1+p_2){|}_k{\mathrm{\Δz}}_{k+1/2}{\mathrm{\Δz}}_k}{E}_x^q{|}_{k+1}\\ =\left(\begin{array}{c}{V}_{Ex}^{q-1}{|}_k+{\mathrm{bV}}_{Ey}^{q-1}{|}_k+\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}({\mathrm{bV}}_{Hx}^{q-1}{|}_{k+1/2}-{\mathrm{bV}}_{Hx}^{q-1}{|}_{k-1/2})\\ +\frac{C_{2z}{|}_k}{(1+p_2){|}_k{\mathrm{\Δz}}_k}(V_{Hy}^{q-1}{|}_{k+1/2}-V_{Hy}^{q-1}{|}_{k-1/2})\end{array}\right)/(1+b^2)\end{array}$

步骤4.2、使用追赶法求解步骤4.1的方程，得到整个计算区域y方向上的电场分量系数

步骤5添加场源到y方向上的磁场分量系数中，并利用步骤3所得x方向上电场分量系数更新计算整个计算区域的y方向上磁场分量系数具体为：

场源的表达式为：

I_my(t)＝exp(-(t-t₀)²/τ²)

式中，t₀,τ为场源参数；

将步骤3所得x方向上电场分量系数代入磁场分量系数计算公式中，更新磁场分量系数

$H_y^q{|}_{k+1/2}=V_{Hy}^{q-1}{|}_{k+1/2}-C_3\frac{C_{2z}{|}_{k+1/2}}{{\mathrm{\Δz}}_{k+1/2}}(E_x^q{|}_{k+1}-E_x^q{|}_k)$

式中

步骤6利用步骤4所得y方向上电场分量系数更新计算整个计算区域的x方向上磁场分量系数具体为

$H_x^q{|}_{k+1/2}=V_{Hx}^{q-1}{|}_{k+1/2}+C_3\frac{C_{2z}{|}_{k+1/2}}{{\mathrm{\Δz}}_{k+1/2}}(E_y^q{|}_{k+1}-E_y^q{|}_k)$

式中

步骤7利用步骤3和步骤4所得电场分量系数，更新计算整个计算区域的极化电流密度具体更新公式为：

$J_x^q{|}_k==0.5{\mathrm{\ϵ}}_{0}s(p_2{E}_x^q{|}_k-p_1{p}_2{E}_y^q{|}_k-p_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_x^k{|}_k+p_1{p}_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_y^k{|}_k)$

$J_y^q{|}_k=0.5{\mathrm{\ϵ}}_{0}s(p_2{E}_y^q{|}_k+p_1{p}_2{E}_x^q{|}_k-p_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_y^k{|}_k-p_1{p}_3\underset{k=0,q>0}{\overset{q-1}{\Σ}}{J}_x^k{|}_k).$

步骤8更新计算整个计算区域的电磁场分量系数的辅助变量，具体更新公式为：

步骤9更新计算观测点处的电磁场分量，具体按照以下公式更新计算：

式中U表示电磁场分量E_x,E_y,H_x,H_y，U^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s＞0的扩展时间，是q阶拉盖尔多项式。

步骤10、将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。

实施例

为了检验本发明方法的正确性和高效性，我们计算了9mm厚碰撞磁化等离子体对垂直入射电磁波的反射和透射系数。入射电磁波为高斯脉冲，加到H_y上，其表达式为exp(-(t-t₀)²/τ²)，式中t₀＝20ps，τ＝5ps。计算区域为350个网格，每个网格大小为75微米，磁化等离子体占据第201到320的网格，其参数为ω_p＝(2π)·50×10⁹rad/s，ω_b＝3.0×10¹¹rad/s，υ＝2.0×10¹⁰Hz，其余为真空。完全匹配层吸收边界放置在计算区域的两端，均为10个网格。时间尺度因子为s＝1.885×10¹²，时间步长为0.25ps,仿真时间为T_f＝1ns，阶数为200。PML吸收边界参数κ_zmax＝1，σ_zmax＝σ_opt，α_zmax＝0。采用本发明方法、WLP-FDFD方法和解析解计算反射系数和透射系数振幅，计算结果如图2和图3。从图中可见，本发明方法和WLP-FDTD方法、解析解计算结果一致，验证了本发明方法的正确性。图4为不同WLP-FDTD方法下以网格数为函数的CPU时间。从图中可以看出因式分解的WLP-FDTD方法计算更为高效。

本发明一种用于磁化等离子体的时域有限差分方法，在直角坐标系下，通过用加权拉盖尔多项式表示电磁场分量，来解时域麦克斯韦方程，使得只是在最后计算观测点处的电磁场分量时用到时间步长，因此计算过程中时间步长可以取得比柯西稳定性条件限制的时间步长更大；在求解电磁场分量系数时，将大型稀疏矩阵方程分裂成两个三对角矩阵方程，使得它在计算时比WLP-FDTD方法更简单、计算速度更快、内存消耗更少而且可以对大区域的电磁场问题进行求解；在设置PML系数时，由于采用了CFS因子，并且通过调整CFS因子中的参数，可以使得该吸收边界对低频与凋落波的吸收更加有效；由于采用了复扩展坐标系，使得PML在实现时避免了场的分裂且与媒质无关。

Claims (10)

1.一种用于磁化等离子体的时域有限差分方法，其特征在于，具体按照以下步骤：

步骤1、输入模型文件；

步骤2、初始化和设置步骤1模型文件中的参数；

步骤3、利用步骤2的参数，更新计算整个计算区域x方向上电场分量系数

步骤4、利用步骤2的参数，更新计算整个计算区域y方向上电场分量系数

步骤5、添加场源到y方向上的磁场分量系数中，并利用步骤3所得x方向上电场分量系数更新计算整个计算区域y方向上的磁场分量系数

步骤6、利用步骤4所得y方向上电场分量系数更新计算整个计算区域的x方向上磁场分量系数

步骤7、利用步骤3和步骤4所得电场分量系数，更新计算整个计算区域的极化电流密度

步骤8、更新计算整个计算区域的电磁场分量系数的辅助变量；

步骤9、更新计算观测点处的电磁场分量；

步骤10、将q+1赋值给q，并判断拉盖尔多项式的阶数q是否达到预设值，若未达到预设值，则返回步骤3，若达到预设值，则结束。

2.根据权利要求1所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤1输入模型文件，具体输入的参数包括：

计算区域大小N_z；空间步长Δz；时间步长Δt；真空中的电导率σ、磁导率μ₀、介电常数ε₀；等离子体中的碰撞频率υ；等离子体频率ω_p；电子回旋频率ω_b；等离子体在计算区域中的位置；吸收边界层数NPML与相关参数κ_zmax，α_zmax，σ_zmax，其中，κ_zmax取整数，κ_zmax取值范围为[1，60]，α_zmax取值范围为[0，1)，σ_zmax/σ_opt取值范围为(0，12]，σ_opt＝(m+1)/150πΔz，m取值范围为[1，20]，Δz取值范围为λ为源的波长；仿真计算时长T_f；加权拉盖尔多项式的阶数q，其中q≥0且为整数；时间尺度因子s，其中s取值范围为[10⁹,10¹³]；观测点；场源参数t₀,τ。

3.根据权利要求2所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤2具体为：

初始化和设置步骤1所得参数：

初始化的参数包括：

将整个计算区域的电磁场分量系数整个计算区域的极化电流密度整个计算区域的电磁场分量系数的和整个计算区域的极化电流密度的和整个计算区域的辅助变量 和拉盖尔多项式 全部初始化为零，其中F_ζ表示E_x,E_y,H_x,H_y；等离子体参数p₁,p₂,p₃初始化为p₁＝0，p₂＝0，PML系数C_1z,C_2z,C₃初始化为C_1z＝1/(1+0.5ε₀s)，C_2z＝1，C₃＝ε₀/μ₀；

设置参数具体为：

设置带有CFS因子的SC-PML吸收边界的参数σ_z,κ_z,α_z：

σ_z＝σ_zmax|z-z₀|^m/d^m

κ_z＝1+(κ_zmax-1)|z-z₀|^m/d^m

α_z＝α_zmax

式中z₀为PML层与非PML截面位置，d是PML吸收边界的厚度；

设置PML系数C_1z,C_2z：

C_1z＝1/(κ_zα_z+σ_z+0.5κ_zε₀s)

C_2z＝(2α_z/ε₀s+1)

设置等离子体参数：

4.根据权利要求3所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤3利用步骤2的参数，更新计算整个计算区域的x方向上电场分量系数具体为：

步骤3.1、将步骤2的参数代入电场分量系数在计算区域的方程，具体为：

式中，k表示第k个计算网格，

步骤3.2、使用追赶法求解步骤3.1的方程，得到整个计算区域x方向上的电场分量系数

5.根据权利要求3所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤4利用步骤2的参数，更新计算整个计算区域的y方向上电场分量系数具体为：

步骤4.1、将步骤2的参数代入电场分量系数在计算区域的方程，具体为：

式中，k表示第k个计算网格，

步骤4.2、使用追赶法求解步骤4.1的方程，得到整个计算区域y方向上的电场分量系数

6.根据权利要求4所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤5添加场源到y方向上的磁场分量系数中，并利用步骤3所得x方向上电场分量系数更新计算整个计算区域的磁场分量系数具体为：

场源的表达式为：

I_my(t)＝exp(-(t-t₀)²/τ²)

式中，t₀,τ为场源参数；

将步骤3所得x方向上电场分量系数代入磁场分量系数计算公式中，更新磁场分量系数

式中

7.根据权利要求5所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤6利用步骤4所得y方向上电场分量系数更新计算整个计算区域的x方向上磁场分量系数具体为

式中

8.根据权利要求1所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤7利用权利要求1中步骤3和步骤4所得电场分量系数，更新计算整个计算区域的极化电流密度具体更新公式为：

9.根据权利要求1所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤8更新计算整个计算区域的电磁场分量系数的辅助变量，具体更新公式为：

10.根据权利要求1所述的一种用于磁化等离子体的时域有限差分方法，其特征在于，所述步骤9更新计算观测点处的电磁场分量，具体按照以下公式更新计算：

式中U表示电磁场分量E_x,E_y,H_x,H_y，U^q表示q阶电磁场分量系数，是q阶加权拉盖尔多项式，是带有时间尺度因子s＞0的扩展时间，是q阶拉盖尔多项式。