Next Article in Journal
Solving Nonlinear Transcendental Equations by Iterative Methods with Conformable Derivatives: A General Approach
Next Article in Special Issue
On Ricci Curvature of a Homogeneous Generalized Matsumoto Finsler Space
Previous Article in Journal
Fractional Equations for the Scaling Limits of Lévy Walks with Position-Dependent Jump Distributions
Previous Article in Special Issue
Harnack Estimation for Nonlinear, Weighted, Heat-Type Equation along Geometric Flow and Applications
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Violation of Neumann Problem’s Solvability Condition for Poisson Equation, Involved in the Non-Linear PDEs, Containing the Reaction-Diffusion-Convection-Type Equation, at Numerical Solution by Direct Method

1
School of Mechanical and Automotive Engineering, South China University of Technology, Guangzhou 510641, China
2
The Faculty of Computational Mathematics and Cybernetics, Lomonosov Moscow State University, 119991 Moscow, Russia
3
Moscow Center of Fundamental and Applied Mathematics, 119234 Moscow, Russia
*
Author to whom correspondence should be addressed.
Submission received: 26 March 2023 / Revised: 18 May 2023 / Accepted: 23 May 2023 / Published: 3 June 2023

Abstract

:
In this paper, we consider the 3D problem of laser-induced semiconductor plasma generation under the action of the optical pulse, which is governed by the set of coupled time-dependent non-linear PDEs involving the Poisson equation with Neumann boundary conditions. The main feature of this problem is the non-linear feedback between the Poisson equation with respect to induced electric field potential and the reaction-diffusion-convection-type equation with respect to free electron concentration and accounting for electron mobility (convection’s term). Herein, we focus on the choice of the numerical method for the Poisson equation solution with inhomogeneous Neumann boundary conditions. Despite the ubiquitous application of such a direct method as the Fast Fourier Transform for solving an elliptic problem in simple spatial domains, we demonstrate that applying a direct method for solving the problem under consideration results in a solution distortion. The reason for the Neumann problem’s solvability condition violation is the computational error’s accumulation. In contrast, applying an iterative method allows us to provide finite-difference scheme conservativeness, asymptotic stability, and high computation accuracy. For the iteration technique, we apply both an implicit alternating direction method and a new three-stage iteration process. The presented computer simulation results confirm the advantages of using iterative methods.

1. Introduction

In the current study, we consider a computer simulation of optical pulse propagation in a semiconductor, which is described by a set of coupled non-linear 3D PDEs. This set involves time-dependent equations describing the evolution of electron and ionized donor concentrations, the equation describing the optical pulse evolution, and the Poisson equation with Neumann boundary conditions (BCs) with respect to the laser-induced electric field potential [1,2,3]. As it is well-known, when solving the Poisson equation, it is necessary to account for the Neumann problem’s solvability condition [4,5], which is a consequence of the Gauss theorem [6]. In the problem under consideration, there is a mutual non-linear influence of the charged particles’ concentrations and the electric field potential because accounting for the electron mobility and describing the evolution of the electron concentration via the reaction-diffusion-convection-type equation. Because of the problem’s nonlinearity, the accuracy of the Poisson equation’s solution is a key issue, and our research aims to address the choices of numerical methods for solving it by comparing the direct method (DM) and the iteration method (IM).
The most popular DM for solving linear elliptic problems stated in the domain of simple (rectangular) geometry is the Fast Fourier Transform (FFT) algorithm. After its development in the second half of the 20th century [7,8], it has been applied ubiquitously in almost every scientific discipline because of its accuracy and cost-effectiveness [9,10,11,12,13,14,15,16,17,18]. Moreover, there are many commercially available packages that implement the FFT algorithm. In the current study, we use Intel MKL Fast Poisson Solver Routines for solving the Poisson equation with Neumann-type inhomogeneous BCs in a 3D case.
However, to our best knowledge, no study has investigated the effectiveness and applicability of DM if the elliptic equation with Neumann BCs is involved in a set of non-linear time-dependent PDEs, which implies the mutual feedback between the Poisson equation solution and a solution of the reaction-diffusion-convection-type equation. As demonstrated below, due to the violation of the Neumann problem’s solvability condition, the DMs (FFT) become ineffective. It is important that the Neumann problem’s solvability condition coincides with the conservation law for the problem under consideration. Therefore, the problem’s invariant is not preserved with high accuracy at using DM, and its value grows over time. Consequently, the Neumann problem’s solvability condition is also violated, and the numerical solution does not converge with those of the differential problem.
To overcome this, we propose to solve the Neumann problem for the Poisson equation using IM. It is known that most IMs are unstable in 3D cases. The stable method for 3D problems, named the implicit alternating direction method (ADM), was proposed almost simultaneously by Peaceman, Rachford, and Douglas [19,20]. It is applied as a finite-difference scheme (FDS) to the parabolic equation and as an iteration technique for solving stationary elliptic problems. Further, in [21], the implicit ADM for 3D problems, which is unconditionally stable and has the second-order approximation both on time coordinates and spatial coordinates, was developed. Let us note that nowadays, ADM and its modifications are extensively used for solving multidimensional non-linear problems [22,23,24,25,26,27,28]. We also propose another IM both for solving the 3D Poisson equation and for the implementation of the conservative non-linear FDS developed for the numerical solution of other equations of the problem based on Crank–Nicolson-type FDS—the original three-stage iteration process (TSIP) [29]. The main advantages of TSIP consist of providing conservativeness property on the iterations and asymptotic stability of FDS. These are crucial features for solving non-linear time-dependent problems.
Below, we compare the effectiveness of using the DM (FFT) and IM (ADM or TSIP) for solving the Neumann problem for the Poisson equation, involved in the set of PDEs, which contains a convection-type term. It should be emphasized that the obtained results and conclusions are regardless of the particular method’s implementation and have a general character because they arise from the solvability condition’s violation caused by the round-off errors accumulation at using DM. It should be emphasized that in [30], the crucial influence of a convection-type term related to electron mobility on the computation’s effectiveness was demonstrated, and at zero-value electron mobility, both methods (DM and IM) possessed the same high accuracy. However, in [30], the problem with homogeneous BCs or a particular case of the external electric field (EEF) action along the direction that is transverse to the optical pulse propagation direction was only considered. On the contrary, below, we study the problem with inhomogeneous BCs corresponding to the action of EEF in one or two directions, including the optical pulse propagation direction.
This paper is organized as follows. In Section 2, we present a mathematical model of the optical pulse propagation in a semiconductor as well as the problem of funding the initial function’s distribution under the condition of the inhomogeneous BCs. Section 3 is devoted to the numerical methods of the problem’s solution. Here, we discuss the DM and IM for a solution to the Neumann problem for the Poisson equation. In Section 4, extensive computer simulation results are given to illustrate the theoretical propositions regarding the effectiveness of the methods for the Poisson equation solution and the influence of EEF on the applicability of these methods. Some concluding remarks are presented in Section 5.

2. Problem Statement

2.1. Semiconductor Plasma Evolution Equations

The 3D problem of semiconductor optical-induced plasma generation is governed by a set of well-known time-dependent, dimensionless, non-linear PDEs [1]:
2 φ x 2 + 2 φ y 2 + 2 φ z 2 = γ ( n N ) , 0 < x < L x ,   0 < y < L y ,   0 < z < L z ,   t > 0 ,
N t = G ( I , δ ) R ( n , N ) , 0 x L x ,   0 y L y ,   0 z L z ,   t > 0 ,
n t = D x x ( n x μ x n φ x ) + D y y ( n y μ y n φ y ) + D z z ( n z μ z n φ z ) + G ( I , δ ) R ( n , N ) , 0 < x < L x ,   0 < y < L y ,   0 < z < L z ,   t > 0 ,
I z + δ 0 δ ( n , N ) I = 0 , 0 x L x ,   0 y L y ,   0 < z L z ,   t > 0 .
Above, the following notations are used: variables x, y, and z are the dimensionless spatial coordinates, and the parameters L x , L y , and L z denote their maximal values, respectively. Variable t represents a dimensionless time unit. The optical pulse propagates along the z-coordinate, and it is named a longitudinal coordinate; the x-coordinate and y-coordinate are transverse ones. Function φ ( x , y , z , t ) is a potential of the optical-induced electric field, function N ( x , y , z , t ) is an ionized donor concentration, function n ( x , y , z , t ) is a free electron concentration in the conduction band of a semiconductor, and function I ( x , y , z , t ) is an intensity of the incident optical pulse. δ ( n , N ) characterizes the coefficient of light energy. Functions G ( I , δ ) and R ( n , N ) describe the generation and recombination of free-charged particles, respectively.
Parameter γ depends on the free-charged particle’s maximal concentration, in particular. Parameter δ 0 denotes a maximal value of the optical energy absorption coefficient. The electron diffusion coefficients D x , D y , and D z and the electron mobility coefficients μ x , μ y , and μ z are positive constants.
A semiconductor is assumed to be a rectangular parallelepiped: Ω = ( 0 < x < L x ;   0 < y < L y ;   0 < z < L z ) bounded by a surface Ω = { 0   ; L x } × { 0   ; L y } × { 0   ; L z } . For the Poisson Equation (1), the Neumann BCs are stated:
φ x | x = 0 , L x = E x , φ y | y = 0 , L y = E y , φ z | z = 0 , L z = E z
The BCs for the Equation (2) are the following:
( n x μ x n φ x ) | x = 0 , L x = ( n y μ y n φ y ) | y = 0 , L y = ( n z μ z n φ z ) | z = 0 , L z = 0
Constants E x , E y , and E z are responsible for the absence (in the case of their zero-value) or presence of EEF along the corresponding spatial coordinate. The conditions (6) correspond to the absence of an electric current through the semiconductor faces.
The intensity distribution of the incident optical radiation has a Gaussian profile of the beam and time-dependent pulse shape:
I | z = 0 = exp ( ( x 0.5 L x 0.1 L x ) 2 ) exp ( ( y 0.5 L y 0.1 L y ) 2 ) ( 1 exp ( 10 t ) ) ,   0 x L x ,   0 y L y , 0 t .
In the initial time moment, the optical pulse intensity is equal to zero:
I ( x , y , z , 0 ) = 0 ,   0 x L x ,   0 y L y ,   0 z L z .
The initial distributions for other functions are defined by the BCs, and they are discussed below in Section 2.2.
The absorption coefficient fundamentally influences the regimes occurring in a semiconductor. Below, we consider its non-linear dependence on free-charged particles concentrations [31,32,33,34]:
δ ( n , N ) = ( 1 N ) e ( ψ ξ n ) ,   ψ = c o n s t > 0 ,   ξ = c o n s t > 0 .
Such an absorption coefficient can result in occurring OB, called a concentration OB. The following functions describe the generation and the recombination of free-charged particles:
G ( I , δ ) = q 0 I δ ( n , N ) ,
R ( n , N ) = n N n e q 2 τ R ,
where a parameter q 0 characterizes the incident optical pulse maximal intensity, the parameter n e q denotes an equilibrium value of the free-charged carrier concentration, and the parameter τ R characterizes a free-electron recombination time.
Remark 1. 
As it is known, the Neumann problem is solved only up to constant, and the multiplicity of its solution occurs. So, changing the electric field potential on a constant does not influence the problem’s (1)–(8) solution if the absorption coefficient is defined by (9). However, if the absorption coefficient depends on the electric field potential, for example, as follows: 
δ ( φ , N ) = ( 1 N ) e α + β | φ | ,   α = c o n s t > 0 ,   β = c o n s t > 0 ,
which takes place if a semiconductor is undergoing the action of a high intensive laser pulse [35,36,37,38], it is necessary to provide re-normalization of the electric field potential computed by solving the Neumann problem for the Poisson equation [39].
We want to stress that due to the presence of electron mobility (coefficients μ x , μ y , μ z ) in the term describing convective transfer in Equation (3), the non-linear feedback between the charged particle concentrations and the electric field potential occurs: the electric field potential influences the free electron concentration evolution and vice versa. Evidently, this non-linear feedback is “switched off” at zero-value electron mobility. In [30], we demonstrated that the presence of such a term plays a crucial role in the Neumann problem’s solvability condition violation using DM.
Let us notice that the existence and uniqueness of the problem (1)–(11) solution was discussed in [40], where the laser pulse propagation was described by the non-linear Schrödinger equation, which is more complicated than Equation (4). We described how to prove the solution’s existence. However, it should be taken into account that the uniqueness of the problem solution can be absent in general cases because the laser beam interaction with the semiconductor under certain conditions possesses the property of optical bistability, which means that for one value of the beam intensity, there are two stable states of electron concentrations and ionized donors’ concentrations.

2.2. The Problem of Finding the Semiconductor Characteristics’ Initial Distributions

For the problem’s numerical solution, it is necessary to compute the initial distributions of the free-electron and ionized-donor concentrations as well as a spatial distribution of the electric field potential. For this purpose, we solve the PDEs (1)–(3) in a steady-state at zero-value optical pulse intensity. In this case, the free-charged carriers’ generation is absent, and the stationary problem is written in the following way:
2 φ i n x 2 + 2 φ i n y 2 + 2 φ i n z 2 = γ ( n i n N i n ) ,
N i n n i n n e q 2 = 0 ,
D x x ( n i n x μ x n i n φ i n x ) + D y y ( n i n y μ y n i n φ i n y ) + D z z ( n i n z μ z n i n φ i n z ) = 0 ,
0 < x < L x , 0 < y < L y , 0 < z < L z ,
with BCs (5)–(6). We denote the solution of the problem (12)–(14), (5)–(6) as:
n i n ( x , y , z ) = n ( x , y , z , t = 0 ) ,   N i n ( x , y , z ) = N ( x , y , z , t = 0 ) ,
φ i n ( x , y , z ) = φ ( x , y , z , t = 0 ) .
Let us note that if the EEF is absent ( E x = E y = E z = 0 ) , then the initial distributions of the functions are homogeneous along the spatial coordinates
φ i n ( x , y , z ) = 0 ,   n i n ( x , y , z ) = N i n ( x , y , z ) = n e q ,   0 x L x ,   0 y L y ,   0 z L z .

2.3. The Charge Conservation Law and Its Correlation with the Neumann Problem’s Solvability

For the problem under consideration, the charge conservation law occurs if the charge flows through the semiconductor surfaces (conditions (6)) are absent.
Q ( t ) = Ω ( n ( x , y , z , t ) N ( x , y , z , t ) ) d x d y d z = c o n s t ,   x , y , z Ω ,   t 0 .
Because a semiconductor is electrically neutral before the optical pulse action, the problem invariant is equal to zero: Q ( t ) = Q ( 0 ) = 0 .
We follow this invariant when providing the computation because it is crucial for a computer simulation of the process: changing the invariant Q ( t ) leads to no correspondence of the numerical solution of the problem to its physical one. In turn, as it is well-known [4,5], the necessary and sufficient condition of the solvability for the Neumann problem for the Poisson equation is the integral equality, which follows from the Gauss theorem [6]:
Ω f ( x , y , z ) d x d y d z = Ω g ( x , y , z ) d x d y d z ,
where f ( x , y , z ) is the right part of the Poisson equation defined in the domain Ω , and g ( x , y , z ) is the right part of Neumann’s BCs. For the problem (1), (5), condition (18) takes the form:
γ Ω ( n ( x , y , z ) N ( x , y , z ) )   d x d y d z = Ω φ v   d x d y d z ,
where
φ v = { φ x ,   x = { 0 , L x } ,   y , z Ω , φ y ,   y = { 0 , L y } ,   x , z Ω , φ z ,   z = { 0 , L z } ,   x , y Ω .
So, the right part of (19) is equal to zero if the BCs (5) are stated, and therefore, the Neumann problem’s solvability condition (18) coincides with the charge conservation law (17). Thus, the FDS’s conservativeness is the necessary condition for providing the Neumann problem (1), (5) solvability condition.

3. Numerical Method

3.1. Conservative Finite-Difference Scheme and Three-Stage Iteration Process for Its Implementation

To solve the set of time-dependent non-linear PDEs in the rectangular 3D domain, we develop the FDS. For this aim, we introduce uniform time and spatial meshes:
ω x = { x i = i h x , i = 0 , P x ¯ , h x = L x / P x } ,   ω y = { y j = j h y , j = 0 , P y ¯ , h y = L y / P y } , ω z = { z k = k h z , k = 0 , P z ¯ , h z = L z / P z } ,   ω z = { z k = ( k 0.5 ) h z , k = 0 , P z + 1 ¯ , h z = L z / P z } , ω t = { t m = m τ , τ = 0 , P t ¯ , τ = L t / P t } ,   ω t = { t m = ( m 0.5 ) τ , m = 1 , P t ¯ , τ = L t / P t } , Ω ^ = ω x × ω y × ω z × ω t ,   Ω ^ = ω x × ω y × ω z × ω t ,
where h x ,   h y ,   h z ,   τ are the mesh steps on the spatial coordinates and on a time coordinate, respectively; P x ,   P y ,   P z ,   P t denote the numbers of mesh nodes on the corresponding coordinates. The mesh functions n h ,   N h ,   φ h are defined on the mesh Ω ^ :
n i j k m = n ( x i , y j , z k , t m ) ,   N i j k m = N ( x i , y j , z k , t m ) ,   φ i j k m = φ ( x i , y j , z k , t m ) .
The mesh function I h is defined on the mesh Ω ^ with shifted nodes both on spatial coordinate z and on time coordinate:
I i j k m = I ( x i , y j , z k , t m ) .
Approximation of the beam intensity on a shifted mesh is caused by the FDS second-order accuracy reaching. For brevity, we use below the index-free notations:
f = f h = f i j k m = f ( x i , y j , z k , t m ) ,   f i ± 1 = f i ± 1 j k m ,   f j ± 1 = f i j ± 1 k m ,   f k ± 1 = f i j k ± 1 m ,   f ^ = f i j k m + 1 ,
where f is one of the mesh functions n h ,   N h ,   φ h , and the following index-free notations for the mesh function I h :
I = I h = I i j k m = I ( x i , y j , z k , t m ) ,   I i ± 1 = I i ± 1 j k m ,   I j ± 1 = I i j ± 1 k m ,   I ^ = I i j k ± 1 m .
We also introduce other notations in the following way:
f i ± 0.5 = 0.5 ( f i + f i + 1 ) , f j ± 0.5 = 0.5 ( f j + f j + 1 ) , f k ± 0.5 = 0.5 ( f k + f k + 1 ) , f 0.5 = 0.5 ( f + f ^ ) = 0.5 ( f i j k m + f i j k m + 1 ) ,   I 0.5 = 0.5 ( I + I ^ ) = 0.5 ( I i j k m + I i j k + 1 m ) , δ = δ h ( n h , N h ) ,   δ ^ = δ ^ h ( n ^ h , N ^ h ) ,   δ 0.5 = 0.5 ( δ + δ ^ ) , G = G h = q 0 I h 0.5 δ ,   G ^ = G ^ h = q 0 I h 0.5 δ ^ ,   G 0.5 = 0.5 ( G + G ^ ) , R = R h = ( n h N h n e q 2 ) / τ R ,   R ^ = R ^ h = ( n ^ h N ^ h n e q 2 ) / τ R ,   R 0.5 = 0.5 ( R + R ^ ) .
The first and the second difference derivatives are defined in a standard way, and they are notated as follows: f x ,   f x ¯ ,   f x ¯ x ,   f y ,   f y ¯ ,   f y ¯ y ,   f z ,   f z ¯ ,   f z ¯ z ,   f t . The 3D Laplace difference operator is also stated in a standard way:
Λ f = f x ¯ x + f y ¯ y + f z ¯ z .
For brevity, we introduce the finite-difference operators:
Ζ x ¯ x ( n , φ ) = μ x h x ( n i + 0.5 φ x n i 0.5 φ x ¯ ) ,   Ζ ^ x ¯ x ( n , φ ) = μ x h x ( n ^ i + 0.5 φ ^ x n ^ i 0.5 φ ^ x ¯ ) , Ζ 0.5 x ¯ x ( n , φ ) = 0.5 ( Z x ¯ x + Z ^ x ¯ x ) ,
and similar operators on other spatial coordinates.
Based on the Crank–Nicolson scheme and using introduced notations, we write the FDS for the problem (1)–(8):
Λ φ ^ = γ ( n ^ N ^ ) ,   i = 1 , P x 1 ¯ , j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ ,   m > 0
N ^ N τ = G 0.5 R 0.5 ,   i = 0 , P x ¯ , j = 0 , P y ¯ ,   k = 0 , P z ¯ ,   m > 0 ,
n ^ n τ = D x ( n 0.5 x ¯ x Ζ 0.5 x ¯ x ( n , φ ) ) + D y ( n 0.5 y ¯ y Ζ 0.5 y ¯ y ( n , φ ) ) + D z ( n 0.5 z ¯ z Ζ 0.5 z ¯ z ( n , φ ) ) + G 0.5 R 0.5 ,
i = 1 , P x 1 ¯ , j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ ,   m > 0 , I ^ I h z + δ 0 δ 0.5 I 0.5 = 0 ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ ,   k = 1 , P z ¯ ,   m > 0 .
The BCs are approximated as follows:
φ ^ x , 0 j k = φ ^ x ¯ , P x j k = E x ,   φ ^ y , i 0 k = φ ^ y ¯ , i P y k = E y ,   φ ^ z , i j 0 = φ ^ z ¯ , i j P z = E z ,
i = 0 , P x ¯ ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ , n ^ x , 0 j k μ x n ^ 0.5 j k φ ^ x , 0 j k = n ^ x ¯ , P x j k μ x n ^ P x 0.5 j k φ ^ x ¯ , P x j k = 0 ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ , n ^ y , i 0 k μ y n ^ i 0.5 k φ ^ y , i 0 k = n ^ y ¯ , i P y k μ y n ^ i P y 0.5 k φ ^ y ¯ , i P y k = 0 ,   i = 0 , P x ¯ ,   k = 0 , P z ¯ , n ^ z , i j 0 μ z n ^ i j 0.5 φ ^ z , i j 0 = n ^ z ¯ , i j P z μ z n ^ i j P z 0.5 φ ^ z ¯ , i j P z = 0 ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ .
The initial intensity distribution is approximated in the following manner:
I i j 0 m + 1 = exp ( ( h x ( i 0.5 P x ) 0.1 L x ) 2 ) exp ( ( h y ( j 0.5 P y ) 0.1 L y ) 2 ) ( 1 exp ( 10 t m + 1 ) )
Let us note that the FDS (20)–(24) has the second order of approximation in the inner mesh nodes, while the BCs are approximated with the first order. This inconsistency is caused by reaching the FDS conservativeness [29]. The difference analog of the invariant (17) is computed by using the formula:
Q ( t m ) = k = 1 P z 1 j = 1 P y 1 i = 1 P x 1 h z h y h x ( n i j k m N i j k m ) .
which possesses the second order of approximation with respect to the mesh node ( x i , y j , z k , t m + 0.5 τ ) , assuming the solution of the differential problem is smooth enough.

3.2. Three-Stage IM for FDS Implementation

As it is well-known, the splitting methods are widely used for the numerical solution of multidimensional linear and non-linear problems. Most of them are effective for 1D and 2D problems, but their generalization for the 3D case is not obvious and often leads to the method’s instability. In [29], we compared the efficiency of the three-stage iterative process and the splitting method and demonstrated the drawbacks of the split-step method for the problem under consideration. The main drawback is a violation of the FDS’s conservativeness. Another one is an absence of asymptotic stability (long-time stability) due to the accumulation of computational error. This manifests in the numerical solution’s symmetry violation contrary to the stated conditions of the process. To overcome this issue, it is necessary to decrease the mesh steps. As the non-stationary process is under consideration, the machine time could become enormous, especially for the 3D problem solution. Contrary to the split-step method, TSIP eliminates these disadvantages, as it provides conservativeness and asymptotic (long-term) stability [29]. In turn, both methods are economic ones.
Below, at the writing iteration stages, we use notations:
Z 0.5 s x ¯ x = 0.5 ( Z x ¯ x ( n , φ ) + Z x ¯ x ( n ^ s , φ ^ s ) ) ,   Z 0.5 s + 1 x ¯ x = 0.5 ( Z x ¯ x ( n , φ ) + Z x ¯ x ( n ^ s + 1 , φ ^ s + 1 ) ) ,   Z 0.5 s + 2 x ¯ x = 0.5 ( Z x ¯ x ( n , φ ) + Z x ¯ x ( n ^ s + 2 , φ ^ s + 2 ) ) , n 0.5 s + 1 x ¯ x = 0.5 ( n x ¯ x + n ^ s + 1 x ¯ x ) ,   I 0.5 s + 1 = 0.5 ( I + I ^ s + 1 ) ,
where s denotes an iteration number. Similar notations are used for the operators Z 0.5 y ¯ y ( n , φ ) ,   Z 0.5 z ¯ z ( n , φ ) as well as for the free-electron concentration and the pulse intensity at the corresponding iteration stages. The stages of the TSIP are the following:
The first stage:
N ^ s + 1 N τ = G 0.5 s , s + 1 R 0.5 s , s + 1 ,
n ^ s + 1 n τ = D x ( n 0.5 s + 1 x ¯ x Ζ 0.5 s x ¯ x ( n , φ ) ) + D y ( n 0.5 s y ¯ y Ζ 0.5 s y ¯ y ( n , φ ) ) + D z ( n 0.5 s z ¯ z Ζ 0.5 s z ¯ z ( n , φ ) ) + G 0.5 s , s + 1 R 0.5 s , s + 1 ,
Λ φ ^ s + 1 = γ ( n ^ s + 1 N ^ s + 1 ) ,
I ^ s + 1 I h z + δ 0 δ 0.5 s + 1 I 0.5 s + 1 = 0 ,
δ 0.5 s + 1 = 0.5 ( δ ( n , N , φ ) + δ ^ ( n ^ s , N ^ s + 1 , φ ^ s ) ) , G ^ s , s + 1 = q 0 I 0.5 s δ ^ ( n ^ s , N ^ s + 1 , φ ^ s ) , R ^ s , s + 1 = ( n ^ s N ^ s + 1 n e q 2 ) / τ R , φ ^ s + 1 x , 0 j k = φ ^ s + 1 x ¯ , P x j k = E x ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ,
( n ^ s + 1 x μ x n ^ s i + 0.5 φ ^ s x , 0 j k ) | i = 0 = ( n ^ s + 1 x ¯ μ x n ^ s i 0.5 φ ^ s x ¯ , P x j k ) | i = P x = 0 ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ , φ ^ s + 1 y , i 0 k = φ ^ s + 1 y ¯ , i P y k = E y ,   i = 0 , P x ¯ ,   k = 0 , P z ¯ ,
( n ^ s + 1 y μ y n ^ s + 1 j + 0.5 φ ^ s + 1 y , i 0 k ) | j = 0 = ( n ^ s + 1 y ¯ μ y n ^ s + 1 j 0.5 φ ^ s + 1 y ¯ , i P y k ) | j = P y = 0 ,   i = 0 , P x ¯ ,   k = 0 , P z ¯ , φ ^ s + 1 z , i j 0 = φ ^ s + 1 z ¯ , i j P z = E z ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ ,
( n ^ s + 1 z μ z n ^ s + 1 k + 0.5 φ ^ s + 1 z , i j 0 ) | k = 0 = ( n ^ s + 1 z ¯ μ z n ^ s + 1 k 0.5 φ ^ s + 1 z ¯ , i j P z ) | k = P z = 0 ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ .
The second stage:
N ^ s + 2 N τ = G 0.5 s + 1 , s + 2 R 0.5 s + 1 , s + 2 ,
n ^ s + 2 n τ = D x ( n 0.5 s + 1 x ¯ x Ζ 0.5 s x ¯ x ( n , φ ) ) + D y ( n 0.5 s + 2 y ¯ y Ζ 0.5 s + 1 y ¯ y ( n , φ ) ) + D z ( n 0.5 s + 1 z ¯ z Ζ 0.5 s + 1 z ¯ z ( n , φ ) ) + G 0.5 s + 1 , s + 2 R 0.5 s + 1 , s + 2 ,
Λ φ ^ s + 2 = γ ( n ^ s + 2 N ^ s + 2 ) ,
I ^ s + 2 I h z + δ 0 δ 0.5 s + 2 I 0.5 s + 2 = 0 ,
δ 0.5 s + 2 = 0.5 ( δ ( n , N , φ ) + δ ^ ( n ^ s + 1 , N ^ s + 2 , φ ^ s + 1 ) ) , G ^ s + 1 , s + 2 = q 0 I 0.5 s + 1 δ ^ ( n ^ s + 1 , N ^ s + 2 , φ ^ s + 1 ) , R ^ s + 1 , s + 2 = ( n ^ s + 1 N ^ s + 2 n e q 2 ) / τ R , φ ^ s + 2 x , 0 j k = φ ^ s + 2 x ¯ , P x j k = E x ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ,
  ( n ^ s + 2 x μ x n ^ s + 2 i + 0.5 φ ^ s + 2 x , 0 j k ) | i = 0 = ( n ^ s + 2 x ¯ μ x n ^ s + 2 i 0.5 φ ^ s + 2 x ¯ , P x j k ) | i = P x = 0 ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ , φ ^ s + 2 y , i 0 k = φ ^ s + 2 y ¯ , i P y k = E y ,   i = 0 , P x ¯ ,   k = 0 , P z ¯ , ( n ^ s + 2 y μ y n ^ s + 1 j + 0.5 φ ^ s + 1 y , i 0 k ) | j = 0 = ( n ^ s + 2 y ¯ μ y n ^ s + 1 j 0.5 φ ^ s + 1 y ¯ , i P y k ) | j = P y = 0 ,   i = 0 , P x ¯ ,   k = 0 , P z ¯ , φ ^ s + 2 z , i j 0 = φ ^ s + 2 z ¯ , i j P z = E z ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ . ( n ^ s + 2 z μ z n ^ s + 2 k + 0.5 φ ^ s + 2 z , i j 0 ) | k = 0 = ( n ^ s + 2 z ¯ μ z n ^ s + 2 k 0.5 φ ^ s + 2 z ¯ , i j P z ) | k = P z = 0 ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ .
The third stage:
N ^ s + 3 N τ = G 0.5 s + 2 , s + 3 R 0.5 s + 2 , s + 3 ,
n ^ s + 3 n τ = D x ( n 0.5 s + 2 x ¯ x Ζ 0.5 s + 2 x ¯ x ( n , φ ) ) + D y ( n 0.5 s + 2 y ¯ y Ζ 0.5 s + 1 y ¯ y ( n , φ ) ) + D z ( n 0.5 s + 3 z ¯ z Ζ 0.5 s + 2 z ¯ z ( n , φ ) ) + G 0.5 s + 2 , s + 3 R 0.5 s + 2 , s + 3 ,
Λ φ ^ s + 3 = γ ( n ^ s + 3 N ^ s + 3 ) ,
I ^ s + 3 I h z + δ 0 δ 0.5 s + 3 I 0.5 s + 3 = 0 ,
δ 0.5 s + 3 = 0.5 ( δ ( n , N , φ ) + δ ^ ( n ^ s + 2 , N ^ s + 3 , φ ^ s + 2 ) ) , G ^ s + 2 , s + 3 = q 0 I 0.5 s + 2 δ ^ ( n ^ s + 2 , N ^ s + 3 , φ ^ s + 2 ) , R ^ s + 2 , s + 3 = ( n ^ s + 2 N ^ s + 3 n e q 2 ) / τ R , φ ^ s + 3 x , 0 j k = φ ^ s + 3 x ¯ , P x j k = E x ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ,
( n ^ s + 3 x μ x n ^ s + 3 i + 0.5 φ ^ s + 3 x , 0 j k ) | i = 0 = ( n ^ s + 3 x ¯ μ x n ^ s + 3 i 0.5 φ ^ s + 3 x ¯ , P x j k ) | i = P x = 0 ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ , φ ^ s + 3 y , i 0 k = φ ^ s + 3 y ¯ , i P y k = E y ,   i = 0 , P x ¯ ,   k = 0 , P z ¯ , ( n ^ s + 3 y μ y n ^ s + 3 j + 0.5 φ ^ s + 3 y , i 0 k ) | j = 0 = ( n ^ s + 3 y ¯ μ y n ^ s + 3 j 0.5 φ ^ s + 3 y ¯ , i P y k ) | j = P y = 0 ,   i = 0 , P x ¯ ,   k = 0 , P z ¯ , φ ^ s + 3 z , i j 0 = φ ^ s + 3 z ¯ , i j P z = E z ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ , ( n ^ s + 3 z μ z n ^ s + 2 k + 0.5 φ ^ s + 2 z , i j 0 ) | k = 0 = ( n ^ s + 3 z ¯ μ z n ^ s + 2 k 0.5 φ ^ s + 2 z ¯ , i j P z ) | k = P z = 0 ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ ,
I i j 0 m + 1 = exp ( ( h x ( i 0.5 P x ) 0.1 L x ) 2 ) exp ( ( h y ( j 0.5 P y ) 0.1 L y ) 2 ) ( 1 exp ( 10 t m + 1 ) ) .
We use values of the mesh functions at the previous time layer as the initial approach for those at the next time layer:
n ^ s = 0 = n ,   N ^ s = 0 = N ,   φ ^ s = 0 = φ .
Non-linear difference Equations (28), (35), and (40), with respect to the free-electron concentration with the corresponding BCs, are solved using the Thomas algorithm. The key issue is the choice of a numerical method for solving the Poisson equation concerning the electric field potential, and it is discussed below in Section 3.3. In turn, the numerical approach for computing the mesh function’s initial distribution is presented in Section 3.4.
The criterion for the iteration process convergence is based on the relative error estimation:
| n ^ s + 3 n ^ s | ε 1 | n ^ s | + ε 2 ,   | N ^ s + 3 N ^ s | ε 1 | N ^ s | + ε 2 ,   | φ ^ s + 3 φ ^ s | ε 1 | φ ^ s | + ε 2 ,   | I ^ s + 3 I ^ s | ε 1 | I ^ s | + ε 2 ,   ε 1 > 0 ,   ε 2 > 0 .
The computation is continued until all inequalities (46) are satisfied in all mesh nodes.

3.3. Numerical Methods for the Neumann Problem Solution

According to the research’s purpose, we solve the Poisson equation at each of the stages by using the DM (FFT) and the IM (ADM). We use the ADM below because the Poisson equation is linear, and, in this case, the three-stage IM may not be efficient in comparison with ADM if the mesh number is not big enough. Then, we compare the computer simulation results obtained using both methods. As mentioned above, the Neumann problem solvability condition coincides with the charge conservation law for the problem under consideration. If we use DM for the Poisson equation solution, then the difference analogs of the ratio (19) should also be valid at the problem’s numerical solution. Evidently, when providing numerical experiments, the conservation law is preserved in the best case within the accuracy of the problem approximation and for the problem under consideration, this is enough for the growth of the invariant’s value because of non-linear feedback between the equations concerning the electron concentration and the electric field potential due to the electron mobility. Therefore, the Neumann problem solvability condition could not be fulfilled. As a consequence, the numerical solution will not approximate a solution to the differential problem.
Consequently, the numerical solution will not approximate a solution of the differential problem.
If we solve the Poisson equation on each of the stages by applying an IM, then it transforms into the difference “heat-conduction equation” by introducing the artificial time step τ ¯ as an iteration parameter:
Λ φ ^ + 1 τ ¯ ( φ ^ φ ) = γ ( n ^ N ^ ) ,   τ ¯ = c o n s t > 0 .
At each of the time layers, one can treat this equation as the difference analog of the differential Helmholtz equation
Δ φ + η φ = f .
The Equation (48) with the Neumann BCs is well known to be uniquely solvable if the following inequality η > 0 is valid [5]. Because the parameter τ ¯ is always positive, the Equation (47) has a unique solution. Consequently, the solvability condition of the difference Equation (47) with the Neumann BCs does not coincide with the conservation law (28) for the problem under consideration in contrast to the case of applying DM to the Poisson equation solution.
To apply the ADM [21], we introduce the auxiliary mesh functions F and Φ , defined on the spatial mesh ω x × ω y × ω z . Function F corresponds to the electric field potential for the iteration process stages (27)–(44): F { φ ^ s + 1 ;   φ ^ s + 2 ;   φ ^ s + 3 } . Function Φ is the corresponding Poisson equation’s right part Φ ^ = γ { ( n ^ s + 1 N ^ s + 1 ) ;   ( n ^ s + 2 N ^ s + 2 ) ;   ( n ^ s + 3 N ^ s + 3 ) } —on each of the iteration stages. As a result, the iteration process for the Poisson equations on each of the stages is written in the following manner:
F b + 1 / 3 F b τ ¯ / 3 = 3 2 F x ¯ x b + 1 / 3 1 2 F b x ¯ x + F b y ¯ y + F b z ¯ z Φ ^ ,
F x b + 1 / 3 | i = 0 = F x ¯ b + 1 / 3 | i = P x = E x , j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ ,
F b + 2 / 3 2 F b + 1 / 3 + F b τ ¯ / 3 = 3 2 ( F y ¯ y b + 2 / 3 F b y ¯ y ) ,
F y b + 2 / 3 | j = 0 = F y ¯ b + 2 / 3 | j = P y = E y , i = 1 , P x 1 ¯ ,   k = 1 , P z 1 ¯ ,
F b + 1 3 / 2 F b + 2 / 3 + 1 / 2 F b τ ¯ / 3 = 3 2 ( F z ¯ z b + 1 F b z ¯ z )
F z b + 1 | k = 0 = F z ¯ b + 1 | k = P z = E z , j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ .
Above b, b + 1 3 , b + 2 3 , b = 0 ,   1 ,   2 are the iteration numbers. As the initial approach of the function F b = 0 , we take the electric field potential computed on the previous iteration of the process (27)–(44).
Excluding F b + 1 / 3 ,   F b + 2 / 3 , the system (49)–(51) could be presented in operator form
F b + 1 F b τ ¯ Λ F b + 1 + F b 2 + τ ¯ 2 4 ( Λ x ¯ x Λ y ¯ y + Λ x ¯ x Λ z ¯ z + Λ y ¯ y Λ z ¯ z ) F b + 1 F b τ ¯ τ ¯ 3 8 Λ x ¯ x Λ y ¯ y Λ z ¯ z F b + 1 F b τ ¯ Φ ^ = 0 .
As it is well-known, such kinds of FDSs are unconditionally stable [21].
Each of the Equations (49)–(51) is solved by using the Thomas algorithm. The process is stopped if the estimation for the residual of the Poisson equation is satisfied:
| F b + 1 x ¯ x + F b + 1 y ¯ y + F b + 1 z ¯ z Φ ^ | < ε 3 ,   ε 3 > 0
for all inner mesh nodes. As a result, the field potential at the corresponding iteration of the multi-stage process is equal to the function F ^ b + 1 (for example, if there is a third stage of the iteration process, then φ ^ s + 3 = F ^ b + 1 ).
Remark 2. 
Alternatively, it is possible to solve the Poisson equation on each of the stages by applying the TSIP:
F ^ b + 1 F τ ¯ = 1 2 ( F ^ b + 1 x ¯ x + F ^ b y ¯ y + F ^ b z ¯ z ) + 1 2 ( F x ¯ x + F y ¯ y + F z ¯ z ) Φ ^ ,
F x b + 1 | i = 0 = F x ¯ b + 1 | i = P x = E x , j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ ,
F ^ b + 2 F τ ¯ = 1 2 ( F ^ b + 1 x ¯ x + F ^ b + 2 y ¯ y + F ^ b z ¯ z ) + 1 2 ( F x ¯ x + F y ¯ y + F z ¯ z ) Φ ^ ,
F y b + 2 | j = 0 = F y ¯ b + 2 | j = P y = E y , i = 1 , P x 1 ¯ ,   k = 1 , P z 1 ¯ ,
F ^ b + 3 F τ ¯ = 1 2 ( F ^ b + 1 x ¯ x + F ^ b + 2 y ¯ y + F ^ b + 3 z ¯ z ) + 1 2 ( F x ¯ x + F y ¯ y + F z ¯ z ) Φ ^ ,
F z b + 3 | k = 0 = F z ¯ b + 3 | k = P z = E z , j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ ,
with the convergence criterion:
| F b + 3 x ¯ x + F b + 3 y ¯ y + F b + 3 z ¯ z Φ ^ | < ε 3 ,   ε 3 > 0 .
Section 4 presents some relevant computer simulation results obtained using TSIP (54)–(57). These computations require more computational time than modified ADM (49)–(51); they have the same number of iterations but different operation numbers.

3.4. Numerical Method for Finding the Initial Distributions of the Mesh Functions

Accounting for the solvability condition of the Neumann problem for the difference Poisson equation, we replace the set of Equations (12)–(14) with equivalent ones by introducing a pseudo-time in each of the equations:
φ t ¯ = 2 φ x 2 + 2 φ y 2 + 2 φ z 2 γ ( n N ) , 0 < x < L x ,   0 < y < L y ,   0 < z < L z ,   t ¯ > 0 ,
N t ¯ = N n n e q 2 r ,
n t ¯ = D x x ( n x μ x n φ x ) + D y y ( n y μ y n φ y ) + D z z ( n z μ z n φ z ) N n n e q 2 r ,
where t ¯ is a positive parameter, which characterizes the artificial time, and constant r is a regularization parameter. The t ¯ -dependent BCs are stated for the electric field potential:
φ ( t ¯ ) x | x = 0 , L x = E x ( 1 e t ¯ τ φ ) , 0 y L y ,   0 z L z ,
φ ( t ¯ ) y | y = 0 , L y = E y ( 1 e t ¯ τ φ ) , 0 x L x ,   0 z L z , φ ( t ¯ ) z | z = 0 , L z = E z ( 1 e t ¯ τ φ ) , 0 x L x ,   0 y L y .
The multiplier ( 1 e t ¯ τ φ ) is introduced to provide smooth reaching a steady-state of the electric field strength, and τ φ is a positive parameter characterizing a rate of exponent tending to zero. The BCs for the free-electron concentration are stated as in (6). The initial conditions for the Equations (58)–(60) are:
φ | t = 0 = 0 ,   n | t = 0 = N | t = 0 = n e q ,   0 x L x ,   0 y L y ,   0 z L z .
So, to obtain initial distributions of the semiconductor characteristics at the arbitrary BCs, we approximate the system (58)–(61), (6) by FDS and realize it by using the iteration process:
N s + 1 N s τ ¯ s = n s N s n 0 2 r ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ,
n s + 1 n s τ ¯ s = D x ( n s + 1 x ¯ x Ζ x ¯ x ( n s , φ s ) ) + D y ( n s y ¯ y Ζ y ¯ y ( n s , φ s ) ) + D z ( n s z ¯ z Ζ z ¯ z ( n s , φ s ) ) n s N s n 0 2 r , i = 1 , P x 1 ¯ ,   j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ , φ s + 1 φ s τ ¯ s = φ s + 1 x ¯ x + φ s y ¯ y + φ s z ¯ z γ ( n s + 1 N s + 1 ) ,   i = 1 , P x 1 ¯ ,   j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ , φ ^ s + 1 x , 0 j k = φ ^ s + 1 x ¯ , P x j k = E x ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ,   ( n ^ s + 1 x μ x n ^ s i + 0.5 φ ^ s x , 0 j k ) | i = 0 = ( n ^ s + 1 x ¯ μ x n ^ s i 0.5 φ ^ s x ¯ , P x j k ) | i = P x = 0 ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ;
N s + 2 N s + 1 τ ¯ s = n s + 1 N s + 2 n 0 2 r ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ,
n s + 2 n s + 1 τ ¯ s = D x ( n s + 1 x ¯ x Ζ x ¯ x ( n s , φ s ) ) + D y ( n s + 2 y ¯ y Ζ y ¯ y ( n s + 1 , φ s + 1 ) ) + D z ( n s + 1 z ¯ z Ζ z ¯ z ( n s + 1 , φ s + 1 ) ) n s + 1 N s + 2 n 0 2 r , i = 1 , P x 1 ¯ ,   j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ , φ s + 2 φ s + 1 τ ¯ s = φ s + 1 x ¯ x + φ s + 2 y ¯ y + φ s + 1 z ¯ z γ ( n s + 2 N s + 2 ) ,   i = 1 , P x 1 ¯ ,   j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ , φ ^ s + 2 y , i 0 k = φ ^ s + 2 y ¯ , i P y k = E y , ( n ^ s + 2 y μ y n ^ s + 1 j + 0.5 φ ^ s + 1 y , i 0 k ) | j = 0 = ( n ^ s + 2 y ¯ μ y n ^ s + 1 j 0.5 φ ^ s + 1 y ¯ , i P y k ) | j = P y = 0 ,   i = 0 , P x ¯ ,   k = 0 , P z ¯ ;
N s + 3 N s + 2 τ ¯ s = n s + 2 N s + 3 n 0 2 r ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ,
n s + 3 n s + 2 τ ¯ s = D x ( n s + 2 x ¯ x Ζ x ¯ x ( n s + 2 , φ s + 2 ) ) + D y ( n s + 2 y ¯ y Ζ y ¯ y ( n s + 1 , φ s + 1 ) ) + D z ( n s + 3 z ¯ z Ζ z ¯ z ( n s + 2 , φ s + 2 ) ) n s + 2 N s + 3 n 0 2 r , i = 1 , P x 1 ¯ ,   j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ , φ s + 3 φ s + 2 τ ¯ s = φ s + 2 x ¯ x + φ s + 2 y ¯ y + φ s + 3 z ¯ z γ ( n s + 3 N s + 3 ) ,   i = 1 , P x 1 ¯ ,   j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ , φ ^ s + 3 z , i j 0 = φ ^ s + 3 z ¯ , i j P z = E z , ( n ^ s + 3 z μ z n ^ s + 2 k + 0.5 φ ^ s + 2 z , i j 0 ) | k = 0 = ( n ^ s + 3 z ¯ μ z n ^ s + 2 k 0.5 φ ^ s + 2 z ¯ , i j P z ) | k = P z = 0 ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ .
The iteration process convergence criterion is given by the formulas:
| n s + 3 N s + 3 n 0 2 | < ε ¯ ,   i = 0 , P x ¯ ,   j = 0 , P y ¯ ,   k = 0 , P z ¯ ,
| D x ( n s + 3 x ¯ x Ζ x ¯ x ( n s + 3 , φ s + 3 ) ) + D y ( n s + 3 y ¯ y Ζ y ¯ y ( n s + 3 , φ s + 3 ) ) + D z ( n s + 3 z ¯ z Ζ z ¯ z ( n s + 3 , φ s + 3 ) ) | < D min ε ¯ , D min = min { D x , D y , D z } ,   i = 1 , P x 1 ¯ ,   j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ , | Λ φ s + 3 γ ( n s + 3 N s + 3 ) | < ε ¯ ,   i = 1 , P x 1 ¯ ,   j = 1 , P y 1 ¯ ,   k = 1 , P z 1 ¯ ,
where ε ¯ is a positive constant. If these inequalities are satisfied, the obtained values are stated as φ i n ,   n i n , and N i n for the iteration process (27)–(44).
Remark 3. 
Above, we proposed the method of computing the initial function distributions if the three components of the EEF are non-zero. However, we should stress that we consider the case of EEF with one or two spatial components in a computer simulation. In the particular case, if the EEF acts only along one spatial coordinate, then the initial function distributions are governed by the following set of equations (for definiteness, we consider   E x 0 ):
2 φ i n ( x , y , z ) x 2 = γ ( n e q e μ x φ i n ( x , y , z ) n e q e μ x φ i n ( x , y , z ) ) , n i n ( x , y , z ) = n e q e μ x φ i n ( x , y , z ) ,   N i n ( x , y , z ) = n e q e μ x φ i n ( x , y , z ) .

4. Computer Simulation Results

Below, we compare the computer simulation results obtained using DM (FFT) and IM for solving the Neumann problem for the Poisson equation involved in the set (1)–(4) of the equations and accounting for the action of EEF on a semiconductor at statement of the BCs. We evaluate the method’s efficiency in dependence on the EEF strength and directions of the electric field action. We also follow the invariant (26) to assess the method’s accuracy.
For the DM implementation, the FFT code belonging to the Intel Math Kernel Library (Poisson Library Routines) is applied to solve the Poisson equation. The solver uses the forward and backward FFT. The description of the library application stresses that “For Poisson equation with Neumann BCs the right-hand side of the problem cannot be arbitrary. Poisson Library verifies that the classical solution exists (up to rounding errors). In any case, Poisson Library computes the normal solution, that is, the solution that has the minimal Euclidean norm of residual” [41].
ADM (49)–(51) as well as TSIPs (27)–(45) are realized as the computer codes developed using C++ language.
The algorithm for solving the problem (20)–(25) is the following. Firstly, we compute the functions at the initial time moment based on the Formulas (63)–(65). Further, we use them to start the iteration process (27)–(44). On each of the iterations, we compute the 3D Poisson equation using one of the discussed methods (ADM or FFT). We provide computations for the following cases: EEF acts only along the longitudinal coordinate or along one of the transverse coordinates (x-axis, for example), or acts simultaneously along longitudinal and transverse ones. The problem (1)–(4) parameters have the following values:
E y = 0 ,   γ = 1000 ,   D x = D y = D z = 10 5 ,   μ x = μ y = μ z = 1 ,   δ 0 = 2 ,   n 0 = 0.01 , q 0 = 1.5 ,   τ R = 1 , ψ = 2 . 553 ,   ξ = 3 , L x = L y = L z = 1 .
The mesh steps and the iteration process parameters are taken as follows:
τ = 10 3 ,   h x = h y = 10 2 ,   h z = 10 2   or   h z = 10 3 ;   ε 1 = 10 6 ,   ε 2 = 10 8 ; τ ¯ = 10 3 ,   ε 3 = 10 3 .
The mesh step on a time coordinate has different values in computational experiments, and it is defined in figure captures. To study the influence of the external electric field component E z , acting along the optical pulse propagation direction, we also provide computation with the decreased mesh step h z . The parameters for the computation of the initial distributions of the functions are chosen as:   τ ¯ s = 10 4 ,     r = 10 3 ,   ε ¯ = 10 4 .
For definiteness, below, discussing the computer modeling results, we mark the electric field potential using the index IM if we solve the Poisson equation using ADM (49)–(51) and the index DM and if we apply the FFT for this aim: φ I M ,   n I M and φ D M ,   n D M , respectively.
In Figure 1, we show 2D projections of the electric field potential and free-electron concentration distributions in a moment of time t = 100. It is easy to see the strong influence of the EEF on them. If the EEF is absent and the incident beam profile is defined in the form (7), then the computed functions are symmetrical with respect to the center of the semiconductor’s illuminated area. Preserving the solution symmetry is an important feature that allows us to estimate the numerical method’s correctness. However, if the EEF acts on a semiconductor along the transverse coordinate, then the function distributions become asymmetric, complicating the evaluation of the obtained numerical solution. The distribution of the free-electron concentration represents a complicated spatial-temporal oscillating structure (Figure 1, right column). It should be stressed that these distributions depend significantly on the choice of numerical method.
As mentioned above, conservativeness is one of the most important properties of the FDS on par with numerical method stability. Let us consider how a Poisson equation’s solution method affects invariant preservation accuracy. As follows from Figure 2, at using IM, the invariant is preserved with very high accuracy: Q ( t ) order is about 10 12 of magnitude if EEF is absent or it acts in one direction (along the optical pulse propagation or in the opposite direction). The invariant value increases by three orders up to a magnitude of 10 8 if the EEF acts in two directions. Thus, the 2D EEF influences the preservation of the invariant more strongly than an action of the 1D EEF. Such Q ( t ) preservation accuracy is enough for providing the conservativeness property of the FDS if accounting for the order of the FDS approximation. Moreover, as follows from Table 1, the invariant value practically does not grow over time after reaching a certain value.
In contrast, at using DM for the Neumann problem solution, the conservativeness of the FDS is violated: the invariant’s value grows over time even at the zero-value EEF (see Figure 3 and Table 1). There is also a dependence of the invariant’s growth on the direction of EEF action. For example, if E x = 0 ,   E z = 10 the invariant is equal to Q ( 100 ) = 1.1 × 10 4 . If the EEF is directed along the z-axis E x = 0 ,   E z = 10 , then the invariant’s value decreases by two orders: Q ( 100 ) = 3.82 × 10 6 . In this case, as follows from Figure 3, the invariant’s growth occurs approximately linearly, and evidently, it continues in time. The reason for such invariant evolution at using DM for the Poisson equation solution is the Neumann problem solvability condition violation. Due to the non-linear feedback between the equations concerning the charged particle concentrations and the electric field potential, the round-off errors accumulate and, consequently, induce errors in computing the right part of the Poisson equation.
Computer experiments show that the EEF presence (a case of the inhomogeneous Neumann BCs) essentially influences the gradient and magnitude of the invariant’s growth. This could be avoided by using a mesh with decreased steps. This can be seen in Figure 3b by comparing the black solid and the red dotted lines. However, decreasing the mesh steps leads to an unacceptable computational cost when solving multidimensional time-dependent problems.
A more detailed quantitative analysis of the invariant change is presented in Table 1. One can see that the invariant’s value at the initial time moment directly depends on the EEF: Q ( 0 ) is equal to zero if the external electric field is absent. If the EEF acts only along the transverse spatial coordinate or only along the longitudinal one, then Q ( 0 ) has an order of magnitude about 10 20 or 10 17 , respectively. The invariant’s value increases up to an order of 10 8 if the EEF acts in both spatial directions. The influence of the choice of the numerical method for the Poisson equation solutions is obvious: we see that invariant Q(t) values with the time growth differ by several orders of magnitude (Table 1). At providing computation using IM, the invariant’s Q ( t ) value increases up to several orders of magnitude (nevertheless maintaining a high accuracy) in the very short time interval at the beginning of computation in comparison with its value at the initial time moment. Then, it practically does not change. The situation becomes the opposite at using DM: regardless of the invariant’s computation accuracy at the initial time moment, its value increases in time.
We stress the absence of the problem solution’s stability if the Poisson equation is solved by using the DM (Table 2). Increasing the computation accuracy for computation of the initial functions’ distribution by increasing the severity of the criterion (66) ( ε ¯ = 10 6 ) does not improve the numerical solution’s stability despite decreasing the invariant by three orders of magnitude at the initial time moment. However, it grows quickly in time and coincides with the Q(t) values obtained with ε ¯ = 10 4 (Table 2, white lines). In contrast, when applying IM, the increase in the accuracy of the initial functions’ distribution computation leads to the corresponding increase in accuracy of the preservation Q(t) (Table 2, gray lines).
Analysis of the numerical solution of the problem (20)–(25) demonstrates a considerable difference between solutions obtained using IM or DM for the Poisson equation solution. Below, we denote them as φ I M ,   φ D M , respectively. Similar notations are used for the free-electron concentrations: n I M ,   n D M . In Figure 4 and Figure 5, we clearly see the difference between the electric field potential φ I M and φ D M distributions along both the x-coordinate and the z-coordinate computed with the same mesh steps. This difference is more significant if the EEF acts in two directions (Figure 5). To understand which solution is more accurate, we provide computations with a refined mesh (Figure 4) and follow the accuracy of invariant Q(t) (Table 3). If the mesh step h z (the step along direction of the EEF acting) is decreased, then the computer simulation results obtained using DM for the Poisson equation solution converge to those obtained by using IM: φ I M and φ D M distributions became much closer (Figure 4, red solid and blue dashed lines).
We see in Table 3 that using DM for the Poisson equation solution with the EEF acting in one direction the decreasing of mesh step h z by one order of magnitude improves the invariant conservation accuracy by two orders of magnitude. However, it is still six orders of magnitude greater than the invariant preservation accuracy reached when using IM. It should be noted that with the same mesh steps, DM computes faster than the IM: computation of up to 100 dimensionless time units takes about 63 h. In contrast, the iteration process requires 190 h. However, decreasing the spatial step in one of the coordinates by an order increases the computation time extremally at using DM: up to 570 h. Computations are implemented on the server Dell Workstation with the processor Intel Xeon 2.6 GHz with 48 cores.
For quantitative assessment of the difference between φ I M and φ D M , we introduce and investigate the following norms:
φ I M φ D M c = max k , j , i | φ I M i j k φ D M i j k | , φ I M φ D M L 2 = ( φ I M φ D M , φ I M φ D M ) = h x h y h z k = 0 P z j = 0 P y i = 0 P x ( φ I M i j k φ D M i j k ) 2 .
Similar norms are used for an assessment of a difference between the free-electron concentrations computed using IM and DM for the numerical solution of the Poisson equation.
As one can see from the graphs depicted in Figure 6, the EEF plays a fundamental role in the DM’s efficiency. If the EEF action is absent, then the C-norm for the difference between φ I M and φ D M has the order of magnitude 10 2 , and the L 2 -norm is about 10 3 (Figure 6a,b). The norm values dramatically grow (up to two orders) if the EEF is non-zero: C-norm value exceeds the unit for some cases. Analyzing Figure 6c,d, we can state that the difference between numerical solutions is maximal in the case of 2D EEF or in the case of EEF acting along longitudinal coordinate. Thus, the Poisson equation solution obtained using DM differs strongly in these cases.
In Figure 7, the free-electron concentration distribution on spatial coordinates is depicted. We see a pronounced difference in both spatial distributions at using IM or DM. They differ both in maximal values of the free-electron concentration and in spatial coordinates of their appearance. We also depict the norms of the difference between functions n I M and n D M values (Figure 8) computed using different numerical methods for the Poisson equation solution. We can see that the distinction between the obtained solutions is significant – it has the order of magnitude 10 1 in the C-norm and 10 3 in the L 2 -norm. Their values depend on the directions of the EEF action and the number of its non-zero components. The L 2 -norm increases exponentially in time (Figure 8b). So, it can be argued that the choice of the computational method for the Poisson equation solution has a fundamental role in the modeling of the process under consideration.
In closing this section, we compare the effectiveness of using two IMs for the Poisson Equation (1) solution: ADM (3.30)–(3.32) and TSIP (54)–(56). To demonstrate the robustness of TSIP, we estimate the computational costs of this process at the computation of the mesh functions at a new time layer. We denote the total number of iterations as S i t , which is defined by summing both the numbers of iteration processes that are necessary to fulfill the conditions (53), (57) and the number of iterations performed to satisfy the criterion (46). As one can see in Figure 9, at the beginning of the computation, the ADM requires significantly fewer iterations than the TSIP. With time increasing, S i t becomes almost the same for both methods. It should be emphasized that the iteration process convergence is strongly determined by the regimes occurring in a semiconductor. Because the physical process under consideration is accompanied by forming complicated contrast dynamical structures, a number of the iterations, which are required for achieving sufficient solution accuracy, can vary significantly according to the problem parameters. From Figure 10, it follows that both methods give numerical solutions that coincide with high accuracy. To assess the difference between the solutions, we introduce additional notations: φ A D M ,   n A D M and φ T S I P ,   n T S I P , which are the functions computed using ADM and TSIP, respectively. The L 1 -norm of the difference between the functions φ A D M ,   φ T S I P is computed as:
φ A D M φ T S I P L 1 = h x h y h z k = 0 P z j = 0 P y i = 0 P x | φ A D M i j k φ T S I P i j k | .
Similar notations are used for the free-electron concentrations. We see that the L 1 -norms of these differences have the order of magnitude 10 7 for the electric field potential (Figure 10a) and 10 9 (Figure 10c) if we provide computations with the iteration parameter τ ¯ = 10 3 and τ ¯ = 10 4 , respectively. The difference between the free-electron concentrations also decreases with the iteration parameter τ ¯ decreasing. Thus, the TSIP method applied for the 3D Poisson equation solution with the Neumann BCs demonstrates robustness and can be used alternatively to the ADM.

5. Conclusions

We demonstrate the inefficiency of DM application for solving the Poisson problem with inhomogeneous Neumann BCs if it is involved to the set of non-stationary non-linear PDEs with realized non-linear feedback. This is a consequence of the violation of the solvability condition of the Neumann problem. We stress that the conservation law of the problem under consideration coincides with the corresponding solvability condition. Therefore, a violation of this law for the FDS leads to a violation of the corresponding solvability condition. Using the IM for the Poisson equation solution allows us to avoid a violation of the Neumann problem solvability condition.
The difference in the results obtained using DM and IM becomes more pronounced at the action of the EEF on a semiconductor. The direction of the EEF action significantly influences the difference between numerical solutions obtained by using IM or DM. When solving the Neumann problem using IM, the problem’s invariant is preserved with high accuracy: it does not grow in time. By decreasing the spatial mesh steps at using DM for a solution of the Poisson equation, it is possible to reach a coincidence of the spatial distribution of the electric field strength obtained using IM with essentially big values of these mesh steps. However, this coincidence occurs during a limited time interval.
We demonstrated that the TSIP (54)–(56) is a promising alternative to the implicit ADM (49)–(51) for a numerical solution of 3D PDEs contained the terms describing reaction-diffusion-convection (in the case under consideration, the generation and recombination of the free charged particles, and electron diffusion and mobility). It possesses high accuracy and asymptotic stability. TSIP is easy to realize and can be used for parallelization algorithms. We highlighted that it is an economical one, as well as the ADM.
We stress that the approach presented in the study could be directly applied to many other scientific problems described by similar PDE sets, for example, in modeling cell chemotaxis or organic molecular systems [42,43,44,45,46].

Author Contributions

Conceptualization, V.T.; methodology, V.T.; software, V.E.; validation, M.L.; formal analysis, V.T. and M.L.; investigation, M.L., V.E. and Y.Y.; writing—original draft preparation, M.L.; writing—review and editing, V.T., M.L., Y.Y. and Z.Y.; visualization, V.E.; supervision, V.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Smith, R.A. Semiconductors; Cambridge University Press: Cambridge, UK, 1978. [Google Scholar]
  2. Markowich, P.A.; Ringhofer, C.A.; Schmeiser, C. Semiconductor Equations; Springer Science & Business Media: Wien, Austria, 2012. [Google Scholar]
  3. Jóźwikowska, A. Numerical solution of the non-linear Poisson equation for semiconductor devices by application of a diffusion-equation finite difference scheme. J. Appl. Phys. 2008, 104, 063715. [Google Scholar] [CrossRef]
  4. Jeffreys, H. Methods of Mathematical Physics, 3rd ed.; Cambridge Mathematical Library: Cambridge, UK, 1972. [Google Scholar]
  5. Atkinson, K.; Hansen, O.; Chien, D. A spectral method for elliptic equations: The Neumann problem. Adv. Comput. Math. 2011, 34, 295–317. [Google Scholar] [CrossRef] [Green Version]
  6. Jackson, J.D. Classical Electrodynamics, 3rd ed.; John Wiley & Sons. Inc.: New York, NY, USA, 1998. [Google Scholar]
  7. Hockney, R.W. A fast direct solution of Poisson’s equation using Fourier analysis. JACM 1965, 12, 95–113. [Google Scholar] [CrossRef]
  8. Cooley, J.W.; Tukey, J.W. An algorithm for the machine calculation of complex Fourier series. Math. Comput. 1965, 19, 297–301. [Google Scholar] [CrossRef]
  9. Dorr, F.W. The direct solution of the discrete Poisson equation on a rectangle. SIAM Rev. 1970, 12, 248–263. [Google Scholar] [CrossRef]
  10. Temperton, C. Direct methods for the solution of the discrete Poisson equation: Some comparisons. J. Comput. Phys. 1979, 31, 1–20. [Google Scholar] [CrossRef]
  11. Buzbe, B.L.; Golub, G.H.; Nielson, C.W. On direct methods for solving Poisson’s equations. SIAM J. Numer. Anal. 1970, 7, 627–656. [Google Scholar] [CrossRef]
  12. Wilhelmson, R.B.; Ericksen, J.H. Direct solutions for Poisson’s equation in three dimensions. J. Comput. Phys. 1977, 25, 319–331. [Google Scholar] [CrossRef]
  13. Schumann, U.; Sweet, R.A. Fast Fourier transforms for direct solution of Poisson’s equation with staggered boundary conditions. J. Comput. Phys. 1988, 75, 123–137. [Google Scholar] [CrossRef] [Green Version]
  14. Boyd, J.P. Chebyshev and Fourier Spectral Methods, 2nd ed.; DOVER Publications, Inc.: New York, NY, USA, 2001. [Google Scholar]
  15. Lai, M.C.; Li, Z.; Lin, X. Fast solvers for 3D Poisson equations involving interfaces in a finite or the infinite domain. J. Comput. Appl. Math. 2006, 191, 106–125. [Google Scholar] [CrossRef] [Green Version]
  16. Wang, T.J.; Sun, T. A spectral method for solving nonhomogeneous Neumann boundary value problems on quadrilaterals. Appl. Numer. Math. 2020, 157, 1–18. [Google Scholar] [CrossRef]
  17. Bergland, G.D. A guided tour of the fast Fourier transform. IEEE Spectr. 1969, 6, 41–52. [Google Scholar] [CrossRef]
  18. Shiferaw, A.; Mittal, R.C. An efficient direct method to solve the three dimensional Poisson’s equation. Am. J. Comput. Math. 2011, 1, 285–293. [Google Scholar] [CrossRef] [Green Version]
  19. Douglas, J. On the Numerical Integration of ∂^2u∂x^2+∂^2u∂y^2=∂u∂t by Implicit Methods. J. Soc. Ind. Appl. Math. 1955, 3, 42–65. [Google Scholar]
  20. Peaceman, D.W.; Rachford, H.H. The numerical solution of parabolic and elliptic differential equations. J. Soc. Ind. Appl. Math. 1955, 3, 28–41. [Google Scholar] [CrossRef]
  21. Douglas, J. Alternating direction methods for three space variables. Numer. Math. 1962, 4, 41–63. [Google Scholar] [CrossRef]
  22. Birkhoff, G.; Varga, R.S. Implicit alternating direction methods. Trans. Am. Math. Soc. 1959, 92, 13–24. [Google Scholar] [CrossRef]
  23. Wachspress, E.L.; Habetler, G.J. An alternating-direction-implicit iteration technique. J. Soc. Ind. Appl. Math. 1960, 8, 403–423. [Google Scholar] [CrossRef]
  24. Marchuk, G.I. Splitting and alternating direction methods. Handb. Numer. Anal. 1990, 1, 197–462. [Google Scholar] [CrossRef]
  25. Zhang, Y.; Sun, Z. Alternating direction implicit schemes for the two-dimensional fractional sub-diffusion equation. J. Comput. Phys. 2011, 230, 8713–8728. [Google Scholar] [CrossRef]
  26. Dai, W.; Nassar, R. Compact ADI method for solving parabolic differential equations. Numer. Methods Part. Differ. Equ. 2002, 18, 129–142. [Google Scholar] [CrossRef]
  27. In’t Hout, K.J.; Welfert, B.D. Stability of ADI schemes applied to convection-diffusion equations with mixed derivative terms. Appl. Numer. Math. 2007, 57, 19–35. [Google Scholar] [CrossRef]
  28. Li, J.; Chen, Y.; Liu, G. High-order compact ADI methods for parabolic equations. Comput. Math. Appl. 2006, 52, 1343–1356. [Google Scholar] [CrossRef] [Green Version]
  29. Trofimov, V.; Loginova, M.; Egorenkov, V. Conservative finite-difference scheme for computer simulation of contrast 3D spatial-temporal structures induced by a laser pulse in a semiconductor. Math. Methods Appl. Sci. 2020, 43, 4895–4917. [Google Scholar] [CrossRef] [Green Version]
  30. Trofimov, V.; Loginova, M.; Egorenkov, V. Numerical methods for solving the 3D Neumann problem of laser-induced plasma evolution in a semiconductor: Direct and iteration methods. Comput. Math. Methods 2020, 3, e1092. [Google Scholar] [CrossRef] [Green Version]
  31. Gibbs, H.M. Optical Bistability: Controlling Light with Light; Academic Press: New York, NY, USA, 1985. [Google Scholar]
  32. Gibbs, H.M.; McCall, S.L.; Venkatesan, T.N.C.; Gossard, A.C.; Passner, A.; Wiegmann, W. Optical bistability in semiconductors. Appl. Phys. Lett. 1979, 35, 451–453. [Google Scholar] [CrossRef]
  33. Smith, S.D. Optical bistability, photonic logic, and optical computation. Appl. Opt. 1986, 25, 1550–1564. [Google Scholar] [CrossRef]
  34. Chai, Z.; Hu, X.; Wang, F.; Niu, X.; Xie, J.; Gong, Q. Ultrafast All-Optical Switching. Adv. Opt. Mater. 2017, 5, 1600665. [Google Scholar] [CrossRef]
  35. Varentsova, S.A.; Trofimov, V.A. The effect of strong light field on the spectrum shift of a Hydrogen-like Atom. Comput. Math. Model. 1999, 10, 275–282. [Google Scholar] [CrossRef]
  36. Varentsova, S.A.; Loginova, M.M.; Trofimov, V.A. Mathematical modeling of optical bistability based on a light-induced electric field. Mosc. Univ. Comput. Math. Cybern. 2003, 1, 22–29. [Google Scholar]
  37. Borovskii, A.V.; Galkin, A.L. Laser Physics: X-ray Lasers, Ultrashort Pulses, Powerful Laser Systems; IzdAT: Moscow, Russia, 1996. [Google Scholar]
  38. Delone, N.B.; Krainov, V.P. Atoms in Strong Light Fields; Springer Series in Chemical Physics; Springer: Berlin/Heidelberg, Germany, 1985. [Google Scholar]
  39. Trofimov, V.A.; Loginova, M.M.; Egorenkov, V.A. A mathematical model of optical bistability and the multiplicity of its solutions. J. Comput. Appl. Math. 2019, 354, 663–681. [Google Scholar] [CrossRef]
  40. Trofimov, V.A.; Loginova, M.M.; Egorenkov, V.A. Conservative finite-difference scheme for the 2D problem of femtosecond laser pulse interaction with kink structure of high absorption in semiconductor. Int. J. Comput. Math. 2020, 97, 207–244. [Google Scholar] [CrossRef]
  41. Available online: https://www.intel.com/ (accessed on 6 March 2019).
  42. Keller, E.F.; Segel, L.A. Model for chemotaxis. J. Theor. Biol. 1971, 30, 225–234. [Google Scholar] [CrossRef] [PubMed]
  43. Murray, J.D. Mathematical Biology II: Saptial Models and Biomedical Application, 3rd ed.; Springer: New York, NY, USA, 2003. [Google Scholar]
  44. Honig, B.; Nicholls, A. Classical electrostatics in biology and chemistry. Science 1995, 268, 1144–1149. [Google Scholar] [CrossRef] [Green Version]
  45. Rodenberger, D.C.; Heflin, J.R.; Garrto, A.F. Excited-state enhancement of optical nonlinearities in linear conjugated molecules. Nature 1992, 359, 309–311. [Google Scholar] [CrossRef]
  46. Li, C.; Deng, X.X.; Wang, Y.X. Non-linear absorption and refraction in multilevel organic molecular system. Chin. Phys. Lett. 2000, 1, 574–576. [Google Scholar] [CrossRef]
Figure 1. Distributions of the electric field potential (a,c,e) and the free-electron concentration (b,d,f) at time moment t = 100 computed using TSIP (27)–(44) for solving the set of Equations (1)–(8) and IM (ADM) (49)–(51) for the Poisson Equation (1) solution with the EEF component values: E x = 0 , E z = 0 (a,c); E x = 10 ,   E z = 10 (b,d); E x = 10 ,   E z = 10 (e,f).
Figure 1. Distributions of the electric field potential (a,c,e) and the free-electron concentration (b,d,f) at time moment t = 100 computed using TSIP (27)–(44) for solving the set of Equations (1)–(8) and IM (ADM) (49)–(51) for the Poisson Equation (1) solution with the EEF component values: E x = 0 , E z = 0 (a,c); E x = 10 ,   E z = 10 (b,d); E x = 10 ,   E z = 10 (e,f).
Mathematics 11 02567 g001
Figure 2. Invariant Q(t) evolution at solving the Poisson Equation (1) using IM (ADM) with the EEF components values: E x = 0 , E z = 0 ((a), red dashed line), E x = 0 , E z = ± 10 ((a), black solid line); E x = 10 , E z = ± 10 ((b), black solid line). The mesh step on z-coordinate h z = 10 2 .
Figure 2. Invariant Q(t) evolution at solving the Poisson Equation (1) using IM (ADM) with the EEF components values: E x = 0 , E z = 0 ((a), red dashed line), E x = 0 , E z = ± 10 ((a), black solid line); E x = 10 , E z = ± 10 ((b), black solid line). The mesh step on z-coordinate h z = 10 2 .
Mathematics 11 02567 g002
Figure 3. Invariant Q(t) evolution at solving the Poisson Equation (1) using DM with the EEF components values: E x = 0 , E z = 0 ,   h z = 10 2 (a, red dotted line), E x = 10 , E z = 0 ,   h z = 10 2 ((a), black solid line), E x = 10 , E z = 10 ,   h z = 10 2 (a, blue dashed line), E x = 0 , E z = 10 ,   h z = 10 2 ((a), blue solid line); E x = 10 , E z = 10 ,   h z = 10 2 ((b), blue dashed line), E x = 0 , E z = 10 , h z = 10 2 ((b), black solid line), E x = 0 , E z = 10 ,   h z = 10 3 ((b), red dotted line).
Figure 3. Invariant Q(t) evolution at solving the Poisson Equation (1) using DM with the EEF components values: E x = 0 , E z = 0 ,   h z = 10 2 (a, red dotted line), E x = 10 , E z = 0 ,   h z = 10 2 ((a), black solid line), E x = 10 , E z = 10 ,   h z = 10 2 (a, blue dashed line), E x = 0 , E z = 10 ,   h z = 10 2 ((a), blue solid line); E x = 10 , E z = 10 ,   h z = 10 2 ((b), blue dashed line), E x = 0 , E z = 10 , h z = 10 2 ((b), black solid line), E x = 0 , E z = 10 ,   h z = 10 3 ((b), red dotted line).
Mathematics 11 02567 g003
Figure 4. The electric field potential distribution along x-coordinate on section y = 0.5 ,   z = 0 (a) and along z-coordinate on section x = 0.5 ,   y = 0.5 (b) in time moment t = 100 computed with the EEF components values E x = 0 ,   E z = 10 using DM with mesh step h z = 10 2 (black dotted line), with mesh step h z = 10 3 (blue dashed line), and IM with mesh step h z = 10 2 (red solid line). Vertical gray dotted lines define the illuminated area of a semiconductor.
Figure 4. The electric field potential distribution along x-coordinate on section y = 0.5 ,   z = 0 (a) and along z-coordinate on section x = 0.5 ,   y = 0.5 (b) in time moment t = 100 computed with the EEF components values E x = 0 ,   E z = 10 using DM with mesh step h z = 10 2 (black dotted line), with mesh step h z = 10 3 (blue dashed line), and IM with mesh step h z = 10 2 (red solid line). Vertical gray dotted lines define the illuminated area of a semiconductor.
Mathematics 11 02567 g004
Figure 5. The electric field potential distribution along x-coordinate in section y = 0.5 ,   z = 0 (a) and along z-coordinate in section x = 0.5 ,   y = 0.5 (b) in a moment t = 100 computed with the EEF components values E x = 10 ,   E z = 10 using IM (red solid line) and DM (black dashed line). The mesh step on z-coordinate h z = 10 2 . Vertical gray dotted lines define the illuminated area of a semiconductor.
Figure 5. The electric field potential distribution along x-coordinate in section y = 0.5 ,   z = 0 (a) and along z-coordinate in section x = 0.5 ,   y = 0.5 (b) in a moment t = 100 computed with the EEF components values E x = 10 ,   E z = 10 using IM (red solid line) and DM (black dashed line). The mesh step on z-coordinate h z = 10 2 . Vertical gray dotted lines define the illuminated area of a semiconductor.
Mathematics 11 02567 g005aMathematics 11 02567 g005b
Figure 6. Evolution of the norms ( C —(a,c); L 2 —(b,d)) of the difference between the electric field potential computed using IM and DM with the EEF components values: E x = 0 , E z = 0 ((a,b)–black solid line), E x = 10 , E z = 0 ((c,d)—red dotted line), E x = 0 , E z = 10 ((c,d)—black solid line); E x = 10 , E z = 10 ((c,d)—blue dashed line). The mesh step on z-coordinate is equal to h z = 10 2 .
Figure 6. Evolution of the norms ( C —(a,c); L 2 —(b,d)) of the difference between the electric field potential computed using IM and DM with the EEF components values: E x = 0 , E z = 0 ((a,b)–black solid line), E x = 10 , E z = 0 ((c,d)—red dotted line), E x = 0 , E z = 10 ((c,d)—black solid line); E x = 10 , E z = 10 ((c,d)—blue dashed line). The mesh step on z-coordinate is equal to h z = 10 2 .
Mathematics 11 02567 g006
Figure 7. The free-electron concentration distribution computed in a moment of time t = 100 along x-coordinate in section y = 0.39 ,   z = 0.09 (a); y = 0.45 ,   z = 0.32 (c) and along z-coordinate in section y = 0.45 ,   y = 0.32 (b); x = 0.61 ,   y = 0.45 with the EEF components values E x = 0 ,   E z = 10 (a,b) and E x = 10 ,   E z = 10 (c,d) at using IM (red solid line) and DM (black dashed line) for the Poisson equation solution. The mesh step on z-coordinate h z = 10 2 .
Figure 7. The free-electron concentration distribution computed in a moment of time t = 100 along x-coordinate in section y = 0.39 ,   z = 0.09 (a); y = 0.45 ,   z = 0.32 (c) and along z-coordinate in section y = 0.45 ,   y = 0.32 (b); x = 0.61 ,   y = 0.45 with the EEF components values E x = 0 ,   E z = 10 (a,b) and E x = 10 ,   E z = 10 (c,d) at using IM (red solid line) and DM (black dashed line) for the Poisson equation solution. The mesh step on z-coordinate h z = 10 2 .
Mathematics 11 02567 g007
Figure 8. Evolution of the norms ( C —(a), L 2 —(b)) of the difference between the free-electron concentrations computed using IM and DM for the EEF components values: E x = 0 ,   E z = 10 ((a,b)—red solid line); E x = 10 ,   E z = 0 ((a,b)—black dotted line); E x = 10 ,   E z = 10 ((a,b)—blue solid line). The mesh step on z-coordinate h z = 10 2 .
Figure 8. Evolution of the norms ( C —(a), L 2 —(b)) of the difference between the free-electron concentrations computed using IM and DM for the EEF components values: E x = 0 ,   E z = 10 ((a,b)—red solid line); E x = 10 ,   E z = 0 ((a,b)—black dotted line); E x = 10 ,   E z = 10 ((a,b)—blue solid line). The mesh step on z-coordinate h z = 10 2 .
Mathematics 11 02567 g008
Figure 9. Number of iterations at solving the Poisson equation using ADM (49)—(51) (black solid line) and TSIP (54)–(56) (red dotted line) and computing with τ ¯ = 10 3 (a) and τ ¯ = 10 4 (b).
Figure 9. Number of iterations at solving the Poisson equation using ADM (49)—(51) (black solid line) and TSIP (54)–(56) (red dotted line) and computing with τ ¯ = 10 3 (a) and τ ¯ = 10 4 (b).
Mathematics 11 02567 g009
Figure 10. L1-norm of difference between solutions obtained using ADM (49)–(51) and TSIP (54)–(56) for the Poisson equation at computing with iteration parameter τ ¯ = 10 3 (a,b) and τ ¯ = 10 4 (c,d): the electric field potential (a,c) and the free-electron concentration (b,d).
Figure 10. L1-norm of difference between solutions obtained using ADM (49)–(51) and TSIP (54)–(56) for the Poisson equation at computing with iteration parameter τ ¯ = 10 3 (a,b) and τ ¯ = 10 4 (c,d): the electric field potential (a,c) and the free-electron concentration (b,d).
Mathematics 11 02567 g010
Table 1. Comparison of the invariant’s Q ( t ) computation accuracy at using IM (ADM) or DM (FFT) for the Poisson equation solution with various values of the external electric field components.
Table 1. Comparison of the invariant’s Q ( t ) computation accuracy at using IM (ADM) or DM (FFT) for the Poisson equation solution with various values of the external electric field components.
E x = 0 E z = 0 E x = 10 E z = 0 E x = 0 E z = 10 E x = 0 E z = 10 E x = 10 E z = 10 E x = 10 E z = 10
t = 009.08 × 1023.73 × 10−173.73 × 10−174.23 × 10−84.23 × 10−8
t = 5IM6.30 × 10−152.68 × 10−122.6810−122.68 × 10−124.23 × 10−84.23 × 10−8
DM9.40 × 10−91.10 × 10−62.27 × 10−61.17 × 10−63.94 × 10−62.78 × 10−6
t = 10IM1.64 × 10−142.69 × 10−122.69 × 10−122.69 × 10−124.23 × 10−84.23 × 10−8
DM5.19 × 10−82.16 × 10−66.32 × 10−61.45 × 10−69.68 × 10−64.77 × 10−6
t = 20IM3.77 × 10−142.71 × 10−122.71 × 10−122.72 × 10−124.23 × 10−84.23 × 10−8
DM3.31 × 10−74.09 × 10−61.67 × 10−57.83 × 10−72.34 × 10−57.47 × 10−6
t = 40IM8.22 × 10−142.76 × 10−122.75 × 10−122.77 × 10−124.23 × 10−84.23 × 10−8
DM8.61 × 10−77.97 × 10−63.95 × 10−52.27 × 10−85.31 × 10−51.34 × 10−5
t = 60IM1.28 × 10−132.81 × 10−122.80 × 10−122.82 × 10−124.23 × 10−84.23 × 10−8
DM1.35 × 10−61.18 × 10−56.32 × 10−55.06 × 10−78.36 × 10−52.05 × 10−5
t = 80IM1.75 × 10−132.86 × 10−122.85 × 10−122.87 × 10−124.23 × 10−84.23 × 10−8
DM1.86 × 10−61.57 × 10−58.74 × 10−51.66 × 10−61.10 × 10−42.82 × 10−5
t = 100IM2.22 × 10−132.91 × 10−122.89 × 10−122.92 × 10−124.23 × 10−84.23 × 10−8
DM2.36 × 10−61.95 × 10−51.10 × 10−43.82 × 10−61.40 × 10−43.66 × 10−5
Table 2. Influence of the accuracy of the initial functions’ distribution on the problem invariant’s computation accuracy computed for the 2D external electric field E x = 10 ,   E z = 10 .
Table 2. Influence of the accuracy of the initial functions’ distribution on the problem invariant’s computation accuracy computed for the 2D external electric field E x = 10 ,   E z = 10 .
t = 0t = 10t = 20t = 60t = 80t = 100
ε ¯ < 10 4 IM4.23 × 10−84.23 × 10−84.23 × 10−84.23 × 10−84.23 × 10−84.23 × 10−8
DM4.23 × 10−89.68 × 10−62.34 × 10−58.36 × 10−51.10 × 10−41.14 × 10−4
ε ¯ < 10 6 IM4.38 × 10−114.38 × 10−114.38 × 10−114.39 × 10−114.40 × 10−114.40 × 10−11
DM4.38 × 10−119.72 × 10−62.35 × 10−58.36 × 10−51.10 × 10−41.14 × 10−4
Table 3. Comparison of invariant Q(t) accuracy at using IM or DM for solving Poisson equation solution with the external electric field E x = 0 ,   E z = 10 and different mesh step h z .
Table 3. Comparison of invariant Q(t) accuracy at using IM or DM for solving Poisson equation solution with the external electric field E x = 0 ,   E z = 10 and different mesh step h z .
DM :   h z = 10 2 DM :   h z = 10 3 IM :   h z = 10 2
t = 52.27 × 10−61.76 × 10−72.68 × 10−12
t = 106.32 × 10−64.8 × 10−72.69 × 10−12
t = 201.67 × 10−51.33 × 10−62.72 × 10−12
t = 403.95 × 10−53.27 × 10−62.77 × 10−12
t = 606.32 × 10−53.27 × 10−62.82 × 10−12
t = 808.74 × 10−57.43 × 10−62.87 × 10−12
t = 1001.10 × 10−49.54 × 10−62.92 × 10−12
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Trofimov, V.; Loginova, M.; Egorenkov, V.; Yang, Y.; Yan, Z. Violation of Neumann Problem’s Solvability Condition for Poisson Equation, Involved in the Non-Linear PDEs, Containing the Reaction-Diffusion-Convection-Type Equation, at Numerical Solution by Direct Method. Mathematics 2023, 11, 2567. https://0-doi-org.brum.beds.ac.uk/10.3390/math11112567

AMA Style

Trofimov V, Loginova M, Egorenkov V, Yang Y, Yan Z. Violation of Neumann Problem’s Solvability Condition for Poisson Equation, Involved in the Non-Linear PDEs, Containing the Reaction-Diffusion-Convection-Type Equation, at Numerical Solution by Direct Method. Mathematics. 2023; 11(11):2567. https://0-doi-org.brum.beds.ac.uk/10.3390/math11112567

Chicago/Turabian Style

Trofimov, Vyacheslav, Maria Loginova, Vladimir Egorenkov, Yongqiang Yang, and Zhongwei Yan. 2023. "Violation of Neumann Problem’s Solvability Condition for Poisson Equation, Involved in the Non-Linear PDEs, Containing the Reaction-Diffusion-Convection-Type Equation, at Numerical Solution by Direct Method" Mathematics 11, no. 11: 2567. https://0-doi-org.brum.beds.ac.uk/10.3390/math11112567

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop