Exponential integrator method for solving the nonlinear Helmholtz equation

1.   Introduction

When the electromagnetic wave propagates in the materials, the medium responses, which reflect the materials' properties such as the magnetic permeability or electric permittivity, usually happen. Considering the propagation of the electromagnetic wave in the nonlinear optics, if one is only interested in the propagation of the monochromatic wave, then the Maxwell's equation which describes the propagation of the electromagnetic wave can be reduced to the nonlinear Helmholtz equation (NLHE) ^{[2,4,10,21,23]} under some reasonable assumptions.

To search efficient numerical methods for solving the NLHE is a very interesting field, especially for the large wave number and the strong nonlinear (large characteristic coefficients in the nonlinear term) problems, which has attracted many attentions in the past decades. There are mainly two topics in the references when designing approximation methods: high accuracy spatial discretization methods and robust iteration schemes for the generated nonlinear equations. For the former, a finite element method was constructed in ^[18] for the problem with discontinuous coefficients. Later, combining a new variable separation method with a fourth order finite difference scheme, an efficient approximation method was investigated in ^[1]. Recently, the existence, uniqueness and the error estimate with explicit wave numbers for the finite element approximation were analyzed in ^[19]. And based on the rearranged Taylor series, we proposed a new finite difference method in ^[10]. As the wave number increases, the solution of the NLHE becomes highly oscillatory, which makes many spatial discratization methods lose their expected accuracy. Besides the approximation method for the spatial discretization, a robust iteration scheme is also necessary to be considered due to the nonlinearity of the NLHE. In ^[2], the Newton's iteration method was investigated. Then, the authors studied the pseudo-time iteration method in ^[21]. Later, a modified Newton's method was also proposed in ^[23]. Recently, we study the error correction iteration method by modifying the original iteration solution with a residual in ^[10]. From the analysis in ^[19], we can see that the convergence of the iteration method heavily depends on the wave number and the characteristic coefficient in the nonlinear term. Moreover, in all methods referred above, linear equations are needed to be solved at each iteration step. However, when the wave number is large, the linear equations generated from the NLHE become indefinite, which are difficult for solving efficiently.

In this paper, we will study the exponential integrator method (EIM) for solving the NLHE. The EIM is a kind of very efficient method for solving the initial-value-type ordinary differential equation, especially for the oscillatory problem ^[7,12]. This idea was extended to solve the second-order oscillatory differential equation in ^[9] and the partial differential equations, such as the nonlinear Dirac equation ^[15,16], the Klein-Gordon equation ^[5,14] and so on. Although the NLHE is a boundary-value-type problem, based on the shooting method, we can transform it into an initial-value-type one. Then the EIM can be applied to the deduced problem, in which the nonlinear term of the NLHE is treated explicitly. Compared with the methods in the references, the proposed method not only can capture the oscillation very well, but also avoid the nonlinear iteration. Moreover, we don't need to solve indefinite linear equations at each iteration step. Therefore, the considered method here can be implemented very efficiently.

The rest of this paper is organized as follows. In the next section, after brief introducing the EIM, we extend the method to approximate the NLHE. Then, in Section 3, we show some numerical examples to confirm the efficiency of the proposed method, including the extension to the nonlinear Helmholtz system. Finally, conclusions are made in Section 4.

2.   Methods for solving the NLHE

In this section, we first briefly introduce the idea of the EIM. Then, based on the shooting method, we will extend it to approximate the NLHE.

2.1.   Exponential integrator method

For completeness, we give a brief introduction to the EIM in this subsection. Considering the second order ordinary differential equation:

where $\Omega$ a positive number, $g$ is a nonlinear function and $p_0, \; q_0$ are two given numbers. We are always interested in the case $\Omega\gg1$. Setting $p(t) = q'(t)$, then (2.1) can be transformed into

Applying the variation-of-constants formula, we arrive at ^[12]

where

Assuming that a uniform partition with a step size $h$ is used for the temporal discretization, i.e., $t_{n} = nh \; (n = 0, 1, 2, \dots),$ approximating $g(q(s))$ with $g(q_n)$ in (2.3), Gautschi ^[8] proposed the one-step method

and the two-step method

where $\text{sinc}\xi = \sin\xi/\xi$, $q_n$ and $p_n$ are the approximations of $q(t_n)$ and $p(t_n)$, respectively.

Then, using the trapezoid rule for $g(q(s))$ in (2.3), Deuflhard ^[6] improved the schemes (2.4) and (2.5) and constructed the folllowing one-step method

and the two-step method

Later, many high-order EIMs are also widely developed, such as the Runge-Kutta EIM ^[20] and the Rosenbrock EIM ^[13,17]. Although we only focus on the simplest schemes (2.4) and (2.5), (2.6) and (2.7) in the following, the idea can be directly extended to these high-order ones.

2.2.   Application to the NLHE

When the electric field $E$ and the material coefficient vary only in one direction $z$, the Maxwell equation can be reduced to the following NLHE under some reasonable assumptions ^{[2,4,10,21,23]}:

where $k_0 = \omega_0/c$ is the wave number with $\omega_0$ being the frequency and $c$ being the speed of light in vacuum, $\varepsilon$ is the nonlinear characteristic coefficient and $i = \sqrt{-1}$ is the imaginary unit. As it is well known that, when the number $k_0$ is large, the Helmholtz equation becomes an indefinite problem and its solution oscillates heavily, which result that the classical numerical methods lose their expected computational effect. The case becomes much worse for the NLHE (2.8) and (2.9) due the existence of the nonlinear effect, which needs to be solved with a robust iteration method.

Next, we apply the EIM introduced above to solve this problem. Rewriting (2.8) and (2.9) as

with $y$ being an initial guess, we transform the NLHE into an initial-value-type problem similar to (2.1) and (2.2) which can be solved by the shooting method to meet the boundary condition at $z = Z_{max}$ in (2.9) in some computational senses. Letting $h = \frac{Z_{max}}{M}$ and $z_m = mh (m = 0, 1, \cdots, M)$, denoting the approximation of $E(z_m)$ by $E_m$ and two different guesses of the shooting method for solving the problem (2.10) by $y_j(j = 0, 1)$ and $p_0^j = (2-y_j)ik_0$, then we can extend the schemes (2.4) and (2.5), (2.6) and (2.7) to solve (2.10) with $\Omega^2 = k_0^2$ and $g = f-k_0^2\varepsilon|E|^2E$. Setting $V^j = p_M^j-ik_0E_M^j$ and the tolerance to be $\delta$, we arrive at the mainly algorithm (a combination of the shooting method and the EIM):

There are mainly three advantages by applying Algorithm 1 to solve the NLHE. First, the oscillatory solution can be captured very well thanks to the computational feature of the EIM. Second, the nonlinear term in the NLHE is explicitly treated in Algorithm 1 which avoids to search a nonlinear iteration method. Third, different from the classical methods to solve indefinite linear equations at each iteration step, only an explicit scheme for an initial-value-type problem need to be solved and lots of computational storages can be saved. Moreover, it is much easier (compared with the finite difference method) to deal with the problem with discontinuous coefficients. Therefore, better computational accuracy and speed are possible by using Algorithm 1.

Remark 1. From the classical theory of the shooting method ^[11], we know that the initial guess greatly impact the numerical result of the shooting method. Generally speaking, an approximation equation's solution (or an approximation of the solution) is a good initial guess when simulating, such as an approximation equation's solution is used in solving the nonlinear Helmholtz equation in the following.

Remark 2. The authors in ^[9,12] have been proved that the EIMs (2.4) and (2.6) are both second-order convergent with respect to the step size $h$. Therefore, if the shooting method is convergent with respect to the initial guesses, Algorithm 1 (SM-EIM) which consists of the shooting method and the EIM will be a second-order scheme, too.

3.   Numerical experiments

In this section, we will show some numerical examples to verify the efficiency of the proposed method above, including the extension to the nonlinear Helmholtz system.

3.1.   Linear equation

First, we consider the following linear Helmholtz equation with a discontinuous coefficient ^[3]:

where $k(z) = \omega/c(z)$ with $\omega$ being a constant and $c(z)$ being a piecewise constant on a partition $0 = z_0 < z_1 < \cdots < z_m = Z_{max}$. Then, this problem's solution is

where the coefficients $\alpha^i_j$ and $\alpha^r_j$ ($\alpha^i_j$ and $\alpha^r_j$ are the imaginary and real parts of $\alpha_j$) are determined by solving a linear system satisfying the $C^1$ compatibility conditions at each point $z_j\; (0 < j < m)$ and the boundary condition $-E'(0) = 1$. For more details, the reader is referred to ^[3].

Letting

and the initial guesses $y_0 = 0$ and $y_1 = 1$ in the shooting method, we test the computational accuracy of Algorithm 1 in $l^\infty$- and $l^2$-norms and collect the results in Tables 1 and 2. Since it is a linear problem (the nonlinear term $g\equiv0$), the SM-EIM has the unique form. As the frequency $\omega$ increasing, the solution of the linear Helmholtz equation becomes more oscillatory, which makes it difficult to be simulated by the classical numerical scheme. But as we can see in Tables 1 and 2, in all tested cases, the SM-EIM achieves very high computational accuracy despite the coefficient is discontinuous, and the results don't change as the frequency developing, which confirm that the proposed method can capture the oscillation solution very well.

3.2.   Nonlinear equation

In this part, we consider the popular example for testing the numerical method's computational effect in the nonlinear Helmholtz equation (see ^[2,10]). For given $f$, $k_0$, $Z_{max}$ and $\varepsilon$, the exact solution of this problem can be determined (see ^[2,10] for the detail). Thus, the errors of the approximation methods can be calculated in this case, too. Setting $f = 0, \; Z_{max} = 10$, the initial guesses $y_0 = 0, \; y_1 = 1$ (two approximation values of the corresponding linear problem) in the shooting method and $\varepsilon = 0.01$, with different $k_0$, we compare the errors of the second-order standard finite difference method (SFDM), the finite volume method (FVM) ^[2], the fourth-order compact finite difference method (CFDM) ^[10] and Algorithm 1 with (2.4) and (2.5) (SM-EIM-G), (2.6) and (2.7) (SM-EIM-D) in Table 3. The results suggest that better simulations can be got by applying the SM-EIM, especially by the high-order approximation method SM-EIM-D. Then, fixed $k_0 = 10$, with other computational parameters taking the same values as above, we show the convergence orders in Figure 1, which is consistent with the claimed second-order convergence in Remark 2. Moreover, fixed $Z_{max} = 1$, we investigate the robustness of different iteration methods, which are performed in Table 4. Compared with that for solving nonlinear equations by the frozen-nonlinearity method (FNM) ^[23], the error correction method (ECM) ^[10], the Newton method (NM) and the modified Newton's method (MNM) ^[23], the robustness for satisfying the shooting accuracy by the SM-EIM-G and SM-EIM-D are much better when the wave number or the characteristic coefficient is large, and this advantage becomes more and more obvious as the parameter increasing.

Finally, we show the numerical solutions got by the SM-EIM-D with $\varepsilon = 0.01, \; 0.1, \; 0.5, \; 1$ and $1.6$ in Figures 2–6, from which we can see that the numerical solutions are almost the same with the exact one. Furthermore, we collect the CPU times(s) in Table 5. It suggests that much less computational time is used in the SM-EIM-D than in the SFDM, and this superiority is enhanced when the step size decreasing. All results verify that the proposed method in this paper is very efficient because of avoiding to search the nonlinear iteration method and to solve indefinite linear equations at each iteration step.

3.3.   Nonlinear system

The third numerical example we consider here is the coupled nonlinear Helmholtz system, which describes the third-harmonic signal generation and enhancement in nonlinear photonic crystals (PhCs) ^[22]

where $z\in(0, Z_{max})$ and the boundary conditions are

where $E_1$ and $E_3$ are the electric field intensities of fundamental frequency and third-harmonic, $z$ is the wave propagation direction, $\alpha_0, \alpha_a, \alpha_b, \gamma_0, \gamma_b$ are constants, $n_i (i = 1, 3)$ is the refractive index, and $\chi_1^{(1)}$ and $\chi_1^{(3)}$ are the third-order nonlinear susceptibility tensors. The system (3.1) and (3.2) can be derived from the Maxwell's equation under the assumption that the electric and magnetic fields are time harmonic and the incident direction of pulsed laser is perpendicular to the 1-D PhCs with a third-order nonlinear medium. For more details, the reader is referred to ^[22]. Besides the difficulty in the NLHE, the difference between the solutions of $E_1$ and $E_3$ is very big due to the third-harmonic intensity $E_3$ are much weaker than the fundamental frequency intensity $E_1$ in this system. Thus, high accurate numerical methods are required when approximating it. Similar to (2.10), rewriting (3.1) and (3.2) to an initial-value-type system and extending Algorithm 1 (SM-EIM) with two nonlinear terms being approximated by $-\frac{3}{4}k_0^2\chi_1^{(3)}\left[\overline{E}_1^{(n)}(E_1^{(n)})^2+E_3^{(n)}(\overline{E}_1^{(n)})^2\right]$ and $-\frac{1}{4}(3k_0)^2\chi_3^{(3)}\left[(E_1^{(n+1)})^3+6E_3^{(n)}E_1^{(n+1)}\overline{E}_1^{(n+1)}\right]$ in sequence, we get the SM-EIM for solving the nonlinear Helmholtz system. Taking $Z_{max} = 1, n_1 = n_3 = \alpha_a = \alpha_b = \alpha_0 = \gamma_b = \gamma_0 = 1, \chi_1^{(3)} = \chi_3^{(3)} = 0.1$, we show the numerical solutions got by applying the SM-EIM (we only used SM-EIM-D here) with $M = 100$ in Figures 7 and 8. And the solutions obtained by the second-order SFDM with $M = 10000$ is used as a reference one. We can find that the numerical solutions generated by the SM-EIM are almost the same with the reference one. When $k_0 = 12$, the imaginary of $E_1$ is about in the range of $-2\times10^{-1}$ to $2\times10^{-1}$, whereas the imaginary of $E_3$ is about in the range of $-2\times10^{-3}$ to $2\times10^{-3}$. These are successfully captured by the proposed method. The similar simulation results are got when $k_0 = 40$. Moreover, we study the CPU time(s) in Table 6, which confirms the fact suggested in the above subsection again. The SM-EIM is also very efficient for solving the nonlinear Helmholtz system.

4.   Conclusions

In this paper, a combination of the shooting method and the exponential integrator method has been investigated for solving the nonlinear Helmholtz equation. The proposed method performs very fast and generates much high accurate simulations. This idea can be extended to other nonlinear problems, too.

Acknowledgments

This work is supported by the Natural Science Foundation of Chongqing (No. cstc2020jcyjmsxmX0551) and Chongqing Key Laboratory of Analytic Mathematics and Applications. We sincerely thank the editor and anonymous referees for their comments that lead to an improvement of this paper.

Conflict of interest

The authors declare no conflict of interest.