CN105760597A - Two-dimensional dispersive medium Crank-Nicolson complete matching layer implementation algorithm based on DG algorithm - Google Patents

Abstract

The invention relates to a two-dimensional dispersive medium Crank-Nicolson complete matching layer implementation algorithm based on a DG algorithm and belongs to the technical field of numerical simulation.The method aims at reducing a two-dimensional dispersive medium FDTD computational domain and simulating limited memory space of a computer into infinite space.The implementation algorithm is technically characterized in that in the process that a two-dimensional modified Maxwell equation with plural stretching coordinate variables is converted into the time domain finite difference from the frequency domain, an auxiliary differential equation method is utilized, and based on the Douglas-Gunn (DG) algorithm, an iteration equation with coefficients being block tridiagonal matrixes is approximately decomposed into two iteration equations with coefficients being tridiagonal matrixes, wherein the two iteration equations can be efficiently solved, and therefore computational efficiency is obviously improved.The implementation algorithm has the advantages of achieving unconditional stability, increasing the electromagnetic field computational speed and saving memory.

Description

Two-dimension chromatic dispersion medium Crank-Nicolson completely permutation based on DG algorithm realizes Algorithm

Technical field

The present invention relates to technical field of value simulation, particularly to one based on Douglas-Gunn (DG) the two-dimension chromatic dispersion medium Crank-Nicolson completely permutation of algorithm realizes algorithm.

Background technology

Time-domain finite difference (FDTD) widely should as a kind of electromagnetic method that calculates In the Electromagnetic Simulation of various time domains calculates, such as antenna, radio circuit, optics and half Conductor etc..FDTD has wide applicability, is suitable for parallel computation, calculation procedure versatility Etc. feature.

But, along with going deep into and various more and more wide variety of demands of scientific research, it is calculated Method originally experience Courant Friedrichs Lewy (CFL) the numerical stability conditions restriction lack Fall into more and more obvious.The suffered numerical stability condition of algorithm itself limits: time during calculating Between step-length and spatial mesh size must be fulfilled for CFL constraints, i.e.

$\mathrm{\Δ}t\≤{\left(c\sqrt{{\left(\mathrm{\Δ}x\right)}^{-2}+{\left(\mathrm{\Δ}y\right)}^{-2}+{\left(\mathrm{\Δ}z\right)}^{-2}}\right)}^{-1}$

In formula, Δ t is for calculating time step, and c is the light velocity of FDTD computational fields medium, Δ x, Δ y and Δ z is three dimensions step-length.In Practical Calculation, spatial spreading step-length and time step Long relative wavelength and cycle are the least, so will necessarily occur when calculating Electrically large size object The situation of inadequate resource, the computational efficiency causing FDTD is the lowest.Therefore to eliminate CFL The restriction of condition, alternating direction implicit (the Alternating Direction of unconditional stability Implicit, ADI) FDTD method, local one-dimension (Local One Dimension, LOD) FDTD method and Crane gram Nicholson (Crank-Nicolson, CN) FDTD method In succession it is suggested.

Although overcoming to a certain extent for ADI-FDTD algorithm and LOD-FDTD algorithm Stability condition restriction, but the computational accuracy of algorithm is too low, and performance is unsatisfactory, its reason It is owing to, after time step increases, the numerical dispersion caused increases, and then causes the mistake of algorithm Difference is bigger.2004, G.Sun et al. used Crank-Nicolson difference scheme to Max Wei Fangcheng carries out sliding-model control, i.e. CN-FDTD, and algorithm is much larger than in time step value Stability condition (such as 20 times) remains to keep good lasting accuracy, shows the suitableeest By property, and CN-FDTD algorithm is a kind of method of easier unconditional stability, will 2 above required in two kinds of algorithms calculating processes are simplified to 1 calculating process, thus significantly Reducing calculation resources, therefore scholars unanimously think that CN-FDTD has broader development Prospect.

Due to the restriction in calculator memory space, numerical computations can only be carried out in limited region, Open or the problem such as the electromagnetic radiation in semi-open region and scattering to be able to simulation, calculating district Absorbing boundary condition is must be provided with, in order to limited mesh space simulation at the cutoff boundary in territory Open infinite space, solves the Electromagnetic Wave Propagation in arbitrary medium and various electromagnetic problem. The completely permutation (Perfectly Matched Layer, PML) proposed by Berenger is mesh The absorbing boundary condition that front application is wider, PML is it is to be understood that pass through in FDTD region A kind of special media layer, the wave impedance of this layer of medium and adjacent media wave resistance are set at cutoff boundary Resist and mate completely, so that incidence wave areflexia ground enters PML layer through interface, PML Layer is lossy dielectric, finally by electro-magnetic wave absorption.The most conventional PML absorbing boundary is main There are stretching coordinate transform completely permutation (SC-PML) and unit anisotropy completely permutation (UPML)。

Owing to the interest of electromagnetic wave with dispersive medium INTERACTION PROBLEMS is continuously increased by people, because of Further investigation is needed in this CN-FDTD and CN-PML emulation relating to dispersive medium badly.

Summary of the invention

It is an object of the invention to for FDTD algorithm by lacking that CFL stability condition is limited Fall into, improve the computational efficiency of PML algorithm of dispersive medium and absorption efficiency and propose based on DG algorithm and Sub-ODE method and the SC-PML algorithm of CN-FDTD.This algorithm The solution efficiency of CN-FDTD-PML algorithm can be significantly improved.

Two-dimension chromatic dispersion Crank-Nicolson completely permutation based on DG algorithm realizes algorithm, Comprise the following steps:

Step 1: Maxwell equation in frequency domain is modified to the Mike with stretching coordinate operator This Wei Fangcheng, and represent in rectangular coordinate system；For dispersive medium item, utilize auxiliary differential Equation method arranges auxiliary variable；

Step 2: according to the mapping transformation relation of frequency domain and time domain, by two in rectangular coordinate system Keep in repair positive Maxwell equation and transform to time-domain representation, be simultaneously based on Sub-ODE method Auxiliary variable is set；

Step 3: time domain expanded form based on Crank-Nicolson Finite Difference Time Domain, Two dimension Maxwell equation in the rectangular coordinate system of forms of time and space is launched into Fdtd Method Form, also auxiliary differential equation is transformed to the form of Fdtd Method simultaneously；

Step 4: Finite Difference-Time Domain form-separating is organized into the form solved, result produces one group The coupling implicit equation in electric field and magnetic field, collated after to obtain coefficient matrix be block tridiagonal matrix The implicit iterative equation of form；

Step 5: use Crank-Nicolson Douglas-Gunn (CNDG) method, will step Electric field iterative equation obtained by rapid 4 be approximately decomposed into can with Efficient Solution, coefficient be three right Two iterative equations of angular moment battle array；

Step 6: utilize the iterative equation obtained by step 5 to solve electric field value, then will solve The electric field value gone out is updated in the iterative equation in magnetic field, solves magnetic-field component, by solve Electric field value and magnetic field value are updated in the iterative equation of auxiliary variable, solve the value of auxiliary variable；

Repetition step 6, thus iterative in time.

Use CNDG method can with effectively by coefficient matrix for block tridiagonal matrix form Electric field iterative equation is approximately decomposed into can with Efficient Solution, coefficient as triple diagonal matrix two Iterative equation, reduces computation complexity, improves computational efficiency, has guidance meaning to FDTD algorithm Justice.

Accompanying drawing illustrates:

Fig. 1 is FB(flow block) of the present invention；

Fig. 2 is the time domain that the present invention blocks Debye dispersive medium under the conditions of different CFLN Relative reflection error characteristics curve；

Fig. 3 is the frequency domain that the present invention blocks Debye dispersive medium under the conditions of different CFLN Relative reflection error characteristics curve；

Fig. 4 is the present invention when blocking Lorentz dispersive medium under the conditions of different CFLN Territory relative reflection error characteristics curve；

Fig. 5 is the frequency that the present invention blocks Lorentz dispersive medium under the conditions of different CFLN Territory relative reflection error characteristics curve.

Detailed description of the invention:

The purport of the present invention is to propose a kind of two-dimension chromatic dispersion medium based on DG algorithm Crank-Nicolson completely permutation realizes algorithm, utilizes Douglas-Gunn to solve thought pole The earth improves Electromagnetic Calculation speed.

Embodiment of the present invention is described further in detail below in conjunction with the accompanying drawings.

Fig. 1 is flow chart of the present invention, implements step as follows:

Step 1: Maxwell equation in frequency domain is modified to the Mike with stretching coordinate operator This Wei Fangcheng, and Maxwell equation revised in frequency domain is represented in rectangular coordinate system； For dispersive medium item, utilize Sub-ODE method that auxiliary variable is set；TM ripple is online Property dispersive medium is propagated and can be described as

$-{\mathrm{j\ωH}}_x=c_0{S}_y^{-1}{\∂}_y{E}_z---\left(1\right)$

${\mathrm{j\ωH}}_y=c_0{S}_x^{-1}{\∂}_x{E}_z---\left(2\right)$

${\mathrm{j\ω\ϵ}}_r\left(\mathrm{\ω}\right)E_z=c_0{S}_x^{-1}{\∂}_x{H}_y-c_0{S}_y^{-1}{\∂}_y{H}_x---\left(3\right)$

In formula, c₀It is free space propagation velocity of electromagnetic wave, S_η(η=x is y) that PML plural number draws Stretch coordinate variable.In the case of PML, S_ηCan be expressed as

$S_{\mathrm{\η}}={\mathrm{\κ}}_{\mathrm{\η}}+\frac{{\mathrm{\σ}}_{\mathrm{\η}}}{{\mathrm{j\ω\ϵ}}_0}---\left(4\right)$

In the case of CFS-PML, S_ηCan be expressed as

$S_{\mathrm{\η}}={\mathrm{\κ}}_{\mathrm{\η}}+\frac{{\mathrm{\σ}}_{\mathrm{\η}}}{{\mathrm{\α}}_{\mathrm{\η}}+{\mathrm{j\ω\ϵ}}_0}---\left(5\right)$

In dispersive medium, ε_r(ω) can be expressed as

${\mathrm{\ϵ}}_r\left(\mathrm{\ω}\right)={\mathrm{\ϵ}}_{\mathrm{\∞}}+\frac{\mathrm{\σ}}{{\mathrm{j\ω\ϵ}}_0}+\mathrm{\χ}\left(\mathrm{\ω}\right)---\left(6\right)$

In formula, ε_∞For radio frequency dielectric constant, σ is electrical conductivity, and χ (ω) is the electric polarization of medium Rate.

Formula (3) can be expressed as

${\mathrm{j\ω\ϵ}}_{\mathrm{\∞}}{E}_z+\frac{\mathrm{\σ}}{{\mathrm{\ϵ}}_0}{E}_z+{\mathrm{j\ωP}}_z=c_0{S}_x^{-1}{\∂}_x{H}_y-c_0{S}_y^{-1}{\∂}_y{H}_x---\left(7\right)$

In formula, P_zCan be obtained by following formula

P_z=χ (ω) E_z (8)

Step 2: according to the mapping transformation relation of frequency domain and time domain, by two in rectangular coordinate system Keep in repair positive Maxwell equation and transform to time-domain representation, be simultaneously based on Sub-ODE method Auxiliary variable is set, i.e.

${\∂}_t{H}_x=-c_0{\mathrm{\κ}}_y^{-1}{\∂}_y{E}_z+g_y---\left(9\right)$

${\∂}_t{H}_y=c_0{\mathrm{\κ}}_x^{-1}{\∂}_x{E}_z-g_x---\left(10\right)$

${\mathrm{\ϵ}}_{\mathrm{\∞}}{\∂}_t{E}_z+\frac{\mathrm{\σ}}{{\mathrm{\ϵ}}_0}\·E_z+{\∂}_t{P}_z=c_0{\mathrm{\κ}}_x^{-1}{\∂}_x{H}_y-c_0{\mathrm{\κ}}_y^{-1}{\∂}_y{H}_x-f_x+f_y---\left(11\right)$

In formula, f_x、f_y、g_xAnd g_yFor auxiliary variable.

Step 3: time domain expanded form based on Crank-Nicolson Finite Difference Time Domain, Two dimension Maxwell equation in the rectangular coordinate system of forms of time and space is launched into Fdtd Method Form, is also transformed to the form of Fdtd Method simultaneously by time domain auxiliary differential equation, utilizes CN item is by formula (9)-(11) discretization, and can obtain discrete equation is

$H_{x_{i,j+1/2}}^{n+1}=H_{x_{i,j+1/2}}^n-C_{{\mathrm{ey}}_{j+1/2}}{\mathrm{\Γ}}_{1y}\left(E_{z_{i,j}}^{n+1}\right)+{\mathrm{\Δt}}_h{\mathrm{\Γ}}_2\left(g_{y_{i,j+1/2}}^{n+1}\right)---\left(12\right)$

$H_{y_{i+1/2,j}}^{n+1}=H_{y_{i+1/2,j}}^n+C_{{\mathrm{ex}}_{i+1/2}}{\mathrm{\Γ}}_{1x}\left(E_{z_{i,j}}^{n+1}\right)-{\mathrm{\Δt}}_h{\mathrm{\Γ}}_2\left(g_{x_{i+1/2,j}}^{n+1}\right)---\left(13\right)$

$\begin{array}{c}{E}_{z_{i,j}}^{n+1}=C_e{E}_{z_{i,j}}^n+C_{{\mathrm{hx}}_i}{\mathrm{\Γ}}_{1x}\left(H_{y_{i+1/2,j}}^{n+1}\right)-C_{{\mathrm{hy}}_j}{\mathrm{\Γ}}_{1y}\left(H_{x_{i,j+1/2}}^{n+1}\right)\\ -\left({\mathrm{\Δt}}_h/{\mathrm{\ϵ}}_{\mathrm{\∞}}\right){\mathrm{\Γ}}_2\left(f_{x_{i,j}}^{n+1}\right)+\left({\mathrm{\Δt}}_h/{\mathrm{\ϵ}}_{\mathrm{\∞}}\right){\mathrm{\Γ}}_2\left(f_{y_{i,j}}^{n+1}\right)-{\mathrm{\ϵ}}_{\mathrm{\∞}}^{-1}{\mathrm{\Γ}}_3\left(P_{z_{i,j}}^{n+1}\right)\end{array}---\left(14\right)$

In formula, ${\mathrm{\Δt}}_h=\frac{\mathrm{\Δ}t}{2},u=2{\mathrm{\ϵ}}_0{\mathrm{\ϵ}}_{\mathrm{\∞}},v=\mathrm{\σ}\mathrm{\Δ}t,C_e=\frac{u-v}{u+v},C_{{\mathrm{e\η}}_{k+1/2}}=\frac{c_0{\mathrm{\Δt}}_h}{{\mathrm{\κ}}_{{\mathrm{\η}}_{k+1/2}}\mathrm{\Δ}\mathrm{\η}},$ $C_{{\mathrm{h\η}}_k}=\frac{{\mathrm{\Δt}}_h}{{\mathrm{\ϵ}}_{\mathrm{\∞}}{\mathrm{\κ}}_{{\mathrm{\η}}_k}\mathrm{\Δ}\mathrm{\η}},r_{1{\mathrm{\η}}_k}=\frac{{\mathrm{\α}}_{{\mathrm{\η}}_k}{\mathrm{\κ}}_{{\mathrm{\η}}_k}+{\mathrm{\σ}}_{{\mathrm{\η}}_k}}{{\mathrm{\κ}}_{{\mathrm{\η}}_k}{\mathrm{\ϵ}}_0}{\mathrm{\Δt}}_h,r_{2{\mathrm{\η}}_k}=\frac{c_0{\mathrm{\Δt}}_h{\mathrm{\σ}}_{{\mathrm{\η}}_k}}{{\mathrm{\ϵ}}_0{\mathrm{\κ}}_{{\mathrm{\η}}_k}^2\mathrm{\Δ}x},L_{1{\mathrm{\η}}_k}=\frac{1-r_{1{\mathrm{\η}}_k}}{1+r_{1{\mathrm{\η}}_k}},$ (η=x, y) is space cell size to Δ η, and (k=i is j) to calculate list to k Insertion numerical value between unit, in order to clear, Г_η[*] is the shorthand in CN method, as

${\mathrm{\Γ}}_{1x}\left(E_{z_{i,j}}^{n+1}\right)=E_{z_{i+1,j}}^{n+1}-E_{z_{i,j}}^{n+1}+E_{z_{i+1,j}}^n-E_{z_{i,j}}^n---\left(15\right)$

For simple linear Debye and Lorentz dispersive medium, χ (ω) can be written as

$\mathrm{\χ}\left(\mathrm{\ω}\right)=\frac{d_0}{e_0+e_{1}j\mathrm{\ω}+e_2{\left(j\mathrm{\ω}\right)}^2}---\left(16\right)$

In formula, d₀、e₀、e₁With e₂For the coefficient of rational polynominal, formula (8) can be discrete for n+1 Formula of time step

$P_{z_{i,j}}^{n+1}=Q_0{P}_{z_{i,j}}^n-Q_1{P}_{z_{i,j}}^{n-1}+Q_e\·\[E_{z_{i,j}}^{n+1}+E_{z_{i,j}}^n+m(E_{z_{i,j}}^n+E_{z_{i,j}}^{n-1})\]---\left(17\right)$

M in formula can be as the switching variable of dispersive medium；Corresponding to Debye and Lorentz look Dispersion media, m takes 0 and 1 respectively；

Step 4: the form of Fdtd Method is organized into the form solved, result produces one Group electric field and the coupled wave equation in magnetic field, this is one group of implicit equation, this prescription journey is decoupled, whole The coefficient electric field implicit iterative equation that the left side is block tridiagonal matrix form is obtained after reason

$(1-D_{2x}-D_{2y})E_{z_{i,j}}^{n+1}=\left(\frac{{\mathrm{\ϵ}}_{\mathrm{\∞}}{C}_e-(1+m)Q_e}{{\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e}+D_{2x}+D_{2y}\right)E_{z_{i,j}}^n+\mathrm{\ρ}---\left(18\right)$

In formula, D_2xIt is defined as

$D_{2x}\left(E_{z_{i,j}}^{n+1}\right)=C_{x_i}^{+}{E}_{z_{i+1,j}}^{n+1}-(C_{x_i}^{+}+C_{x_i}^{-})E_{z_{i,j}}^{n+1}+C_{x_i}^{-}{E}_{z_{i-1,j}}^{n+1}---\left(19\right)$

D_2yThere is similar definition；ρ is field amount and the shorthand of auxiliary variable known to the n moment, It is defined as

$\begin{array}{c}\mathrm{\ρ}={\mathrm{\ξ}}_{x_i}\·\left(H_{y_{i+1/2,j}}^n-H_{y_{i-1/2,j}}^n\right)-{\mathrm{\ζ}}_{y_j}\·\left(H_{x_{i,j+1/2}}^n-H_{x_{t,j-1/2}}^n\right)\\ -S_{x_i}^{+}{g}_{x_{i+1/2,j}}^n+S_{x_i}^{-}{g}_{x_{t-1/2,j}}^n-S_{y_j}^{+}{g}_{y_{i,j+1/2}}^n+S_{y_j}^{-}{g}_{y_{i,j-1/2}}^n-{\mathrm{\θ}}_{x_i}{f}_{x_{i,j}}^n\\ -\frac{Q_0-1}{{\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e}{P}_{z_{i,j}}^n+\frac{Q_1}{{\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e}{P}_{z_{i,j}}^{n-1}-\frac{{\mathrm{mQ}}_e}{{\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e}{E}_{z_{i,j}}^{n-1}+{\mathrm{\θ}}_{y_j}{f}_{y_{t,j}}^n\end{array}---\left(20\right)$

Step 5: formula (18) remains a more complicated matrix, it is still desirable to the biggest Amount of calculation, uses CNDG method, adds respectively on the left side of formula (18) and the rightWithArrange

$(1-D_{2x})E_{z_{i,j}}^{*}=\left(\frac{{\mathrm{\ϵ}}_{\mathrm{\∞}}{C}_e-(1+m)Q_e}{{\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e+D_{2x}+2D_{2y}}\right)E_{z_{i,j}}^n+\mathrm{\ρ}---\left(21\right)$

$(1-D_{2y})E_{z_{i,j}}^{n+1}=E_{z_{i,j}}^{*}-D_{2y}{E}_{z_{i,j}}^n---\left(22\right)$

Step 6: utilize the iterative equation obtained by step 5 to solve electric field value, then will solve The electric field value gone out is updated in the iterative equation in magnetic field, solves magnetic-field component, by solve Electric field value and magnetic field value are updated in the iterative equation of auxiliary variable, solve the value of auxiliary variable.

Repetition step 6, thus iterative in time.

Fig. 2 is the time domain that the present invention blocks Debye dispersive medium under the conditions of different CFLN Relative reflection error characteristics curve, Fig. 4 is that the present invention blocks under the conditions of different CFLN The time domain relative reflection error characteristics curve of Lorentz dispersive medium, the extracting method in order to verify, Inventive algorithm is programmed, obtains result shown in attached 2 and Fig. 4 by Computer Simulation, Wherein,In formulaIt is that conventional FDTD algorithm ensure that numerical value Discrete interval maximum time of stability.As can be seen from Figure, the absorbent properties of CN-PML Do not change with the increase of CFLN；Fig. 3 is that the present invention blocks under the conditions of different CFLN The frequency domain relative reflection error characteristics curve of Debye dispersive medium, Fig. 5 is that the present invention is in difference The frequency domain relative reflection error characteristics curve of Lorentz dispersive medium is blocked under the conditions of CFLN, Can be seen that PML absorbent properties do not change with the increase of CFLN, illustrate that this algorithm has Having unconditional stability, simulation process required time is shorter compared with traditional algorithm simulation time.

The foregoing is only presently preferred embodiments of the present invention, be not limiting as the present invention, all at this Within bright spirit and principle, any modification, equivalent substitution and improvement etc. made, all should wrap Within being contained in protection scope of the present invention.

Claims (4)

1. two-dimension chromatic dispersion medium Crank-Nicolson completely permutation based on DG algorithm realizes Algorithm, comprises the following steps:

Step 1: Maxwell equation in frequency domain is modified to the Maxwell with stretching coordinate operator Equation, and represent in rectangular coordinate system；For dispersive medium item, utilize auxiliary differential equation Method arranges auxiliary variable；

Step 2: according to the mapping transformation relation of frequency domain and time domain, by two maintenances in rectangular coordinate system Positive Maxwell equation transforms to time-domain representation, is simultaneously based on Sub-ODE method and arranges Auxiliary variable；

Step 3: time domain expanded form based on Crank-Nicolson Finite Difference Time Domain, by time In the rectangular coordinate system of territory form, two dimension Maxwell equation is launched into the shape of Fdtd Method Formula, is also transformed to the form of Fdtd Method simultaneously by auxiliary differential equation；

Step 4: Finite Difference-Time Domain form-separating is organized into the form solved, result produces one group of electric field With the coupling implicit equation in magnetic field, collated after to obtain coefficient matrix be block tridiagonal matrix form Implicit iterative equation；

Step 5: use Crank-Nicolson Douglas-Gunn method, by the electricity obtained by step 4 Iterative equation be approximately decomposed into can with Efficient Solution, coefficient as triple diagonal matrix two repeatedly For equation；

Step 6: utilize the iterative equation obtained by step 5 to solve electric field value, then will solve Electric field value is updated in the iterative equation in magnetic field, solves magnetic-field component, the electric field that will solve Value and magnetic field value are updated in the iterative equation of auxiliary variable, solve the value of auxiliary variable；

Repetition step 6, thus iterative in time.

2. complete according to the two-dimension chromatic dispersion medium Crank-Nicolson based on DG algorithm described in right 1 Matching layer realizes algorithm, it is characterised in that: step 2, the Maxwell equation conversion that will revise To time domain, auxiliary variable f is set simultaneously_x、f_y、g_xAnd g_y。

3. according to the two-dimension chromatic dispersion medium Crank-Nicolson based on DG algorithm described in right 1 Completely permutation realizes algorithm, it is characterised in that: step 3, based on Crank-Nicolson time domain The time domain expanded form of finite-difference algorithm

$H_{x_{i,j+1/2}}^{n+1}=H_{x_{i,j+1/2}}^n-C_{{\mathrm{ey}}_{j+1/2}}{\mathrm{\Γ}}_{1y}\left(E_{z_{i,j}}^{n+1}\right)+{\mathrm{\Δt}}_h{\mathrm{\Γ}}_2\left(g_{y_{i,j+1/2}}^{n+1}\right)---\left(1\right)$

$H_{y_{i+1/2,j}}^{n+1}=H_{y_{i+1/2,j}}^n+C_{{\mathrm{ex}}_{i+1/2}}{\mathrm{\Γ}}_{1x}\left(E_{z_{i,j}}^{n+1}\right)-{\mathrm{\Δt}}_h{\mathrm{\Γ}}_2\left(g_{x_{i+1/2,j}}^{n+1}\right)---\left(2\right)$

$\begin{array}{c}{E}_{z_{i,j}}^{n+1}=C_e{E}_{z_{i,j}}^n+C_{{\mathrm{hx}}_i}{\mathrm{\Γ}}_{1x}\left(H_{y_{i+1/2,j}}^{n+1}\right)-C_{{\mathrm{hy}}_j}{\mathrm{\Γ}}_{1y}\left(H_{x_{i,j+1/2}}^{n+1}\right)\\ -({\mathrm{\Δt}}_h/{\mathrm{\ϵ}}_{\mathrm{\∞}}){\mathrm{\Γ}}_2\left(f_{x_{i,j}}^{n+1}\right)+({\mathrm{\Δt}}_h/{\mathrm{\ϵ}}_{\mathrm{\∞}}){\mathrm{\Γ}}_2\left(f_{y_{i,j}}^{n+1}\right)-{\mathrm{\ϵ}}_{\mathrm{\∞}}^{-1}{\mathrm{\Γ}}_3\left(f_{z_{i,j}}^{n+1}\right)\end{array}---\left(3\right)$

In formula, Γ_η[*] is the shorthand in CN method, as

${\mathrm{\Γ}}_{1x}\left(E_{z_{i,j}}^{n+1}\right)=E_{z_{i+1,j}}^{n+1}-E_{z_{i,j}}^{n+1}+E_{z_{i+1,j}}^n-E_{z_{i,j}}^n---\left(4\right)$

4. complete according to the two-dimension chromatic dispersion medium Crank-Nicolson based on DG algorithm described in right 1 Matching layer realizes algorithm, it is characterised in that: step 4, the iterative equation of electric field component is carried out Arrange, can be to obtain the coefficient electric field implicit iterative equation that the left side is block tridiagonal matrix form

$\begin{array}{c}(1-D_{2x}-D_{2y})E_{z_{i,j}}^{n+1}\\ =\{\[{\mathrm{\ϵ}}_{\mathrm{\∞}}{C}_e-(1+m)Q_e\]/({\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e)+D_{2x}+D_{2y}\}{E}_{z_{i,j}}^n+\mathrm{\ρ}\end{array}---\left(5\right)$

In formula, ρ is field amount and the shorthand of auxiliary variable known to the n moment, is defined as

$\begin{array}{c}\mathrm{\ρ}={\mathrm{\ξ}}_{x_i}\·(H_{y_{i+1/2,j}}^n-H_{y_{i-1/2,j}}^n)-{\mathrm{\ζ}}_{y_j}\·(H_{x_{i,j+1/2}}^n-H_{x_{i,j-1/2}}^n)\\ -S_{x_i}^{+}{g}_{x_{i+1/2,j}}^n+S_{x_i}^{-}{g}_{x_{i-1/2,j}}^n-S_{y_j}^{+}{g}_{y_{i,j+1/2}}^n+S_{y_j}^{-}{g}_{y_{i,j-1/2}}^n-{\mathrm{\θ}}_{x_i}{f}_{x_{i,j}}^n\end{array}$

$-\frac{Q_0-1}{{\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e}{P}_{z_{i,j}}^n+\frac{Q_1}{{\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e}{P}_{z_{i,j}}^{n-1}-\frac{{\mathrm{mQ}}_e}{{\mathrm{\ϵ}}_{\mathrm{\∞}}+Q_e}{E}_{z_{i,j}}^{n-1}+{\mathrm{\θ}}_{y_j}{f}_{y_{i,j}}^n---\left(6\right)$