Exponentially fitted multisymplectic scheme for conservative Maxwell equations with oscillary solutions

Aiming at conservative Maxwell equations with periodic oscillatory solutions, we adopt exponentially fitted trapezoidal scheme to approximate the temporal and spatial derivatives. The scheme is a multisymplectic scheme. Under periodic boundary condition, the scheme satisfies two discrete energy conservation laws. The scheme also preserves two discrete divergences. To reduce computation cost, we split the original Maxwell equations into three local one-dimension (LOD) Maxwell equations. Then exponentially fitted trapezoidal scheme, applied to the resulted LOD equations, generates LOD multisymplectic scheme. We prove the unconditional stability and convergence of the LOD multisymplectic scheme. Convergence of numerical dispersion relation is also analyzed. At last, we present two numerical examples with periodic oscillatory solutions to confirm the theoretical analysis. Numerical results indicate that the LOD multisymplectic scheme is efficient, stable and conservative in solving conservative Maxwell equations with oscillatory solutions. In addition, to one-dimension Maxwell equations, we apply least square method and LOD multisymplectic scheme to fit the electric permittivity by using exact solution disturbed with small random errors as measured data. Numerical results of parameter inversion fit well with measured data, which shows that least square method combined with LOD multisymplectic scheme is efficient to estimate the model parameter under small random disturbance.


Introduction
Maxwell equations are basic and important mathematical physical models in electromagnetism. They reflect the general law of electromagnetic field excited by charge current and the internal motion of electromagnetic field. The study of the equations is helpful to deepen and enrich the understanding of the materiality of electromagnetic field. Efficient numerical analysis to Maxwell equation is an important study hotspot. Numerical investigations of Maxwell equations include symplectic exponential integrator [1], compact scheme [2], Runge-Kutta method [3], stochastic Euler scheme [4], multisymplectic schemes [5][6][7][8] method [9], wavelet collocation method [10,11], Hamiltonian time integrators [12], finite element method [13,14], spectral method [15] and other numerical methods.
In the past few decades, multisymplectic schemes, as a kind of structure-preserving algorithms, are popular and gain wide applications in scientific computing due to long-time stability in numerical simulation [5-8, 11, 16-19]. A review on stochastic multisymplectic schemes for Maxwell equations was given by Zhang et al. in [5]. Cui et al. [17] presented stochastic multisymplectic method to solve random Schrödinger equation. Hong et al. [6,7] proposed multisymplectic methods and investigated discrete conservative quantity to stochastic Maxwell equations driven by additive noise. Multisymplectic wavelet collocation method for Maxwell's equations were considered in [10,11]. Zhu et al. [18] also applied multisymplectic wavelet collocation methods to solve Schrödinger equations. High-order compact multisymplectic method was presented to simulate Schrödinger equations in [19]. Kong et al. [8] also studied splitting multisymplectic integrators for Maxwell's equations. Since much attention has been gained in multisymplectic schemes, we are committed to the investigations on multisymplectic schemes to simulate Maxwell equations with oscillatory solutions.
However, the investigation of multisymplectic scheme for oscillatory solutions of Maxwell equations is extremely sparse. In this paper, we apply exponentially fitted trapezoidal scheme to construct multisymplectic scheme to simulate Maxwell equations with oscillatory solutions. Our main contributions are highlighted as follows: • We adopt exponentially fitted method to establish a multisymplectic scheme for solving Maxwell equations with periodic solutions.
• To reduce the calculations, we present a LOD multisymplectic scheme by means of splitting method.
• We rigorously prove the conservation, stability, dispersion and convergence theorems of our method.
• Numerical experiment of Maxwell equations verifies the effectiveness of the LOD multisymplectic scheme in simulating high-frequency oscillatory equations. Meanwhile numerical inversion results show that least square method and LOD multisymplectic scheme can be combined to fit the electric permittivity under small measured errors.
The rest of this paper is outlined as follows. In section 2, we commence by reviewing some preliminary knowledge of Maxwell equations and recall exponentially fitted trapezoidal scheme. In section 3, a new multisymplectic scheme is proposed by utilising exponentially fitted trapezoidal scheme and we proceed with the proof of multisymplectic conservation law and two energy conservation laws. Then we built a splitting multisymplectic scheme by splitting original equations into three LOD systems in section 4. We also analyze its numerical stability, dispersion relation and convergence. In section 5, we present two numerical examples to verify the theoretical analysis and confirm the effectiveness of our method. Numerical test is also given on fitting the electric permittivity by applying the least square method and exponentially fitted trapezoidal scheme.

Preliminaries
Maxwell equations are the theoretical basis of electrodynamics. Maxwell equations are equivalent to multisymplectic Hamiltonian systems. From the multisymplectic form, we can construct multisymplectic schemes of Maxwell equations.

Multisymplectic structure of Maxwell equations
For isotropic homogeneous medium, by assuming that there is no free charge and conduction current, the propagation law of electromagnetic field can be described by Maxwell's equations in the following form of curl Here E = (E x , E y , E z ) T is electric field and H = (H x , H y , H z ) T is magnetic field. ε and μ denote permittivity and permeability, respectively. Below ε and μ are constants.
The corresponding component form of the system (1) is as follows By introducing a new variable Z = (H x , H y , H z , E x , E y , E z ) T , the system (2) is equivalent to the following multisymplectic Hamiltonian system Here the Hamiltonian function is S(Z) = 0, and corresponding antisymmetric matrices are Under periodic boundary condition, Maxwell equations have two energy invariants as follows: In the case of one-dimensional transverse magnetic polarization, the electric field is E = (0, 0, E z ) T , and the magnetic field is H = (0, H y , 0) T . Then the transmission law of electromagnetic field accords with the following one-dimensional Maxwell equations Exponentially fitted trapezoidal scheme The exponentially fitted Runge-Kutta method is a kind of efficient numerical tool to simulate ordinary differential equation with oscillatory solutions [27][28][29][30]. For example, exponentially fitted trapezoidal scheme is a special exponentially fitted Runge-Kutta method. Yin et al. [30] applied exponentially fitted trapezoidal scheme to simulate a kind of stochastic oscillator system and obtained efficient and stable numerical results. For initial value problem of ordinary differential equation u 0 = f(t, u), u(t 0 ) = u 0 , exponentially fitted trapezoidal scheme has the following formula where ω t denotes the frequency parameter of the scheme [31]. By Taylor expansion, the local truncation error of the method is ðu ð3Þ þ o 2 t u 0 ÞDt 3 =12: Therefore, exponentially fitted trapezoidal scheme is convergent with second order accuracy. If we choose ω t such that u ð3Þ þ o 2 t u 0 ¼ 0; third order convergence is obtained.

Exponentially fitted multisymplectic scheme of Maxwell equations
Based on the multisymplectic formula (3), we can design multisymplectic schemes of Maxwell equations with oscillatory solutions by applying exponentially fitted Runge-Kutta method. Below, we take exponentially fitted trapezoidal scheme as an example.

Derivation of exponentially fitted multisymplectic scheme
In temporal and spatial discretization of multisymplectic Hamiltonian system (3), we apply exponentially fitted trapezoidal scheme. First, divide the solution area into equidistant meshes. The step sizes taken in the directions of time and space are denoted by Δt, Δx, Δy, Δz, respectively. Denote the mesh points in the directions of time and space by The numerical solutions of u(t, x, y, z) at points ðt n ; x i ; y j ; z k Þ; ðt n ; x; y; zÞ; ðt; x i ; y; zÞ; ðt; x; y j ; zÞ; ðt; x; y; z k Þ are denoted by u n i;j;k ; u n ; u i ; u j ; u k ; respectively. We also make the following notions: Define the difference operator as follows Below, we apply exponentially fitted trapezoidal scheme in temporal and spatial directions to discretize the multisymplectic Hamiltonian system (3) and get the numerical method as follows: For convenience, we give the following simple notions in above formula: where ω t , ω x , ω y , ω z are the frequency parameters of exponentially schemes applied in the direction of time and space, respectively. For the sake of simplicity, below we omit the semi-node indices. The omitted indices in below are the semi-node indices. We abbreviate above formula as The corresponding component form of above scheme (9) is Numerical properties of exponentially fitted multisymplectic scheme Next we analyze and give numerical properties of our proposed scheme (9). Theorem 1 The scheme (9) satisfies the following multisymplectic conservation law Here for the sake of simplicity, we omit the semi-node indices. So the scheme (9) is a multisymplectic scheme of Maxwell equations.
Proof Taking variation on both sides of the formula (9) yields that Then by making wedge product on both sides of above formula with dZ nþ 1 2 and considering the properties of antisymmetric matrix and wedge product, we obtain the discrete multisymplectic conservation law (11). Below, we refer to the scheme (9) as exponentially fitted multisymplectic (abbreviated to EFMS) scheme.
Theorem 2 Under periodic boundary condition, EFMS scheme (9) has two discerete energy invariants as follows: and other two norms have similar definitions. Here for the sake of simplicity, we omit the semi-node indices.
Proof By taking inner product on both sides of EFMS scheme (9) with Z nþ 1 2 we obtain that Above formula can be rewritten as the component form as follows where the omitted indices are semi-nods indices. Obviously In view of periodic boundary condition, we get that and other five similar equalities. Therefore, summing all terms on both sides of above equation (13) with respect to the indices i, j, k yields the first discrete energy conservation law I n 1 ¼ I nþ1 1 : By applying difference operator δ t to EFMS scheme (9), we obtain that Making inner product on both sides of above formula with d t Z nþ 1 2 yields that where the omitted indices are semi-nods indices. Clearly, On account of periodic boundary condition and commutativity of difference operators, we derive that and other five similar equalities. Hence, by summing all terms on both sides of above equation (14) with respect to the indices i, j, k, we obtain the second discrete energy conservation law I n 2 ¼ I nþ1 2 : As is well known that if the media is lossless, the system (1) is divergence-free [5,6], i.e., The divergence convergence is satisfied by EFMS scheme. Theorem 3 EFMS scheme (9) preserves the following two discrete divergences: Here the discrete divergence operator is defined by and other symbols are similar defined. (16) can be seen as the discrete conservation law of (15). Proof Above operator definition yields that According to formula (10) of EFMS scheme, we obtain that and other two similar equalities, which completes the proof of the first divergence convergence. The proof of the second assertion is similar.

LOD exponentially fitted multisymplectic scheme
According to the splitting method of vector field [32][33][34], to reduce calculation cost, we consider the following LOD system of multisymplectic Hamiltonian system (3) The system (17) has LOD multisymplectic conservation law We recast the system (17) into the following six one-dimensional Maxwell equations: Using exponentially fitted trapezoidal scheme to solve numerically the LOD system (17) yields the following LOD scheme Here for the sake of simplicity, we omit the semi-node indices in above formula. Similarly to Theorem 1, we can prove that above scheme satisfies the discrete LOD multisymplectic conservation law as follows: So we call the scheme (22) LOD exponentially fitted multisymplectic (abbreviated to LODEFMS) scheme. LODEFMS scheme (22) is equivalent to the following six one-dimensional scheme: Theorem 4 LODEFMS scheme (22) is unconditionally stable and non-dissipative scheme.
Proof For simplicity, we consider the following one-dimensional Maxwell equations and corresponding LODEFMS scheme: To analyze the stability of scheme (28), define the following wave solution where � i; o x denote imaginary unit and the wave number, respectively. ½E 0 y ; H 0 z � T is a nonzero eigenvector of scheme (28), and ρ is the stability factor. Inserting (29) into scheme (28) yields that ðr À 1Þ cos y À � i 3a t εa x sin yðr þ 1Þ À � i 3a t ma x sin yðr þ 1Þ ðr À 1Þ cos y (34). Taking the following equivalence relation into consideration yields that, as the step sizes tend to 0, numerical dispersive relations (32)- (34) converge to the exact relations (36), respectively. We can see that in the two cases numerical errors of phase velocity are very small. Similar results in other data cases can be obtained.
Theorem 6 Under precise initial condition, numerical solutions of LODEFMS scheme (22) converge to solutions of Maxwell equations with the following numerical errors: where r and s are error vectors of electric field and magnetic field, respectively. Proof For simplicity, we consider errors of one-dimensional scheme to solve one-dimensional Maxwell equation (27). According to local truncation analysis of scheme (28), we obtain that where r n i ¼ E y ðx i ; t n Þ À ðE y Þ n i ; s n i ¼ H z ðx i ; t n Þ À ðH z Þ n i ; Multiplying the first and second equation of (39) with r ; respectively, yields that ε a t r nþ1 By adding both sides of above equation and summing with respect to the subscript i, we can get Applying Cauchy inequality leads to So we derive by recursion that According to Gronwall inequality, we obtain the following error estimation: Under precise initial condition, (41) reduces to (38). By observing above inequality, we can find that, as the initial data are sufficiently accurate, LODEFMS scheme is convergent.

Numerical examples
In this section, we apply our scheme LODEFMS scheme (22) to two Maxwell equations with periodic boundary condition and oscillatory solutions, to veirfy our theoretical analysis.

One-dimensional example
We consider one-dimensional Maxwell equations (7) with the following oscillatory solution where ε γ = k β, μ β 2 = ε. Without loss of generality, here we take two sets of test data as follows:  Next analyse the conservation results of LODEFMS scheme. Fig 5 plots two discrete energy I n 1 and I n 2 in first testing case. The curves in the figure are similar to horizontal straight line. Error evolution curve at adjacent time of two discrete energy in second testing case is depicted in Fig 6. The error magnitude of I n 1 and I n 2 are approximately 10 −14 and 10 −13 , respectively. This verifies that I n 1 and I n 2 are two conservation quantity of LODEFMS scheme. Least square method is an effective tool to solve inverse problem of parameter estimation [35]. Next, we adopt least square method and LODEFMS scheme to fit the parameter ε in one-dimensional Maxwell equations. Measured data in simulation is theoretical solution with random disturbance δ. Suppose that random perturbation δ obeys a normal distribution with expectation 0 and variance σ. Now given simulation data at points where E k zi ; H k yi are our numerical solutions. Assume that the real parameter in experiment is ε = 1.   Table 1 lists the inverting parameter ε � under different distribution σ. Figs 7 and 8 show measured data and inversion data with x = π, t 2 [0, 10] and x 2 [0, 2π], t = 0.5, respectively. The results show that the inversion data is consistent with measured data. The inverting parameters estimated by least square method and LODEFMS scheme are reasonable under small disturbance error.

PLOS ONE
From the figures, we observe that the errors of LODEFMS scheme are smaller than those of LOD central box scheme. Next we check energy conservation property of LODEFMS scheme. Generally, we take step sizes as Dt ¼ 0:04; Dx ¼ Dy ¼ Dz ¼ p 4k : Figs 13 and 14 show errors of discrete energy I n 1 and I n 2 in two testing cases, respectively. We find that, the error magnitudes of I n 1 and I n 2 are about 10 −15 and 10 −12 , respectively. Same conservation results can be obtained by other step sizes. This verifies that LODEFMS scheme preserves two discrete energy conservation laws.
At last, to test the conservation of discrete divergence, we define the errors of two discrete divergences by Err-DivEðnÞ ¼ DxDyDz X i;j;k j � r i;j;k � E nþ1 À � r i;j;k � E n j; Err-DivHðnÞ ¼ DxDyDz X i;j;k j � r i;j;k � H nþ1 À � r i;j;k � H n j;  In other data cases, we can obtain the same conservation behavior of discrete divergences.

Conclusion
Conservative Maxwell equations with periodic oscillatory solutions are investigated numerically. First we have constructed a conservative multisymplectic scheme, preserving two discrete energy invariants and two discrete divergences, by applying exponentially fitted trapezoidal scheme to solve Maxwell equations with periodic solutions. Then LOD multisymplectic scheme is presented in order to reduce the cost of calculation. Unconditional stability, dispersion and convergence analysis for LODEFMS scheme are established additionally. We have carried out two examples in simulating Maxwell equations with periodic oscillatory solutions to illustrate the effectiveness and conservation property of LODEFMS scheme. In addition, taking one-dimensional Maxwell equation as an example, we find that least square method and LODEFMS scheme can be combined to fit the electric permittivity under small random perturbation data.