Conservative finite‐difference scheme for computer simulation of contrast 3D spatial‐temporal structures induced by a laser pulse in a semiconductor

The numerical approach for computer simulation of femtosecond laser pulse interaction with a semiconductor is considered under the formation of 3D contrast time‐dependent spatiotemporal structures. The problem is governed by the set of nonlinear partial differential equations describing a semiconductor characteristic evolution and a laser pulse propagation. One of the equations is a Poisson equation concerning electric field potential with Neumann boundary conditions that requires fulfillment of the well‐known condition for Neumann problem solvability. The Poisson equation right part depends on free‐charged particle concentrations that are governed by nonlinear equations. Therefore, the charge conservation law plays a key role for a finite‐difference scheme construction as well as for solvability of the Neumann difference problem. In this connection, the iteration methods for the Poisson equation solution become preferable than using direct methods like the fast Fourier transform. We demonstrate the following: if the finite‐difference scheme does not possess the conservatism property, then the problem solvability could be broken, and the numerical solution does not correspond to the differential problem solution. It should be stressed that for providing the computation in a long‐time interval, it is crucial to use a numerical method that possessing asymptotic stability property. In this regard, we develop an effective numerical approach—the three‐stage iteration process. It has the same economic computing expenses as a widely used split‐step method, but, in contrast to the split‐step method, our method possesses conservatism and asymptotic stability properties. Computer simulation results are presented.


| INTRODUCTION
Investigation of a high-intensity laser pulse interaction with a semiconductor is a very modern problem, and many authors referring to it. [1][2][3][4][5][6][7][8] Semiconductors are widely used due to the various nonlinear responses under a laser radiation action. Under certain conditions, the semiconductor response becomes bistable. This phenomenon allows the hysteresis loop implementing for the semiconductor characteristics representing existence of two stable states (an upper state and a low one) for a semiconductor response at the same incident optical pulse intensity. Therefore, the possibility of switching between these states at different intensities can be realized, and all-optical data processing based on the optical bistability (OB) phenomenon is developed. [9][10][11][12][13][14][15][16] From a mathematical point of view, the OB problem is described by a set of nonlinear time-dependent partial differential equations (PDEs). As a rule, it is necessary to solve this set of PDEs numerically using a computer simulation because of the complicated spatiotemporal structures formation with strong gradients on their boundaries (contrast structures) under a laser pulse action. Then, special attention should be paid to the accuracy of the difference solution. Another characteristic feature of this problem arises from the presence of the Poisson equation among this PDEs set. As a rule, the Neumann boundary conditions (BCs) or third-order BCs accompanies the Poisson equation. As it is wellknown, for the Neumann problem solvability, the BCs and the right part of the Poisson equation must satisfy the definite conditions, [17][18][19][20] which are necessary to take into account when we solve the Poisson equation numerically. For solving a single Poisson equation in a simple domain, the direct methods, [21][22][23] for example, fast Fourier transform (FFT), are usually used. The application of this method possesses well-known advantages. Much rare the iteration method is used for this aim. So, in Trofimov and Filipchuk, 24 we solve the Poisson equation in a complicated domain using the iteration process. In the case under consideration below, the right part of the Poisson equation depends on the solution of the other two nonlinear equations describing the evolution of the charged particle concentrations. Thus, because of the medium nonlinear feedback, the computational errors arising from solving one of the PDEs lead to computational errors accumulation and distortion of the right part of other PDEs. By turn, this results in the violation of the Neumann problem solvability condition. Because this condition depends on the conservation law for the equations set, then the finite-difference scheme (FDS) conservatism property plays a key role in computer simulation of the problem under consideration.
In our research, we compare two approaches for the Poisson equation solution with Neumann BCs: the direct method based on FFT and the iteration process on the base of the split-step scheme (stabilization method). The paper computer simulation results demonstrate that the last one is much preferable and that the direct method has an essential disadvantage.
As it is well-known, splitting methods are widely used for solving multi-dimensional nonlinear problems. [25][26][27][28][29][30][31][32] Therefore, in Trofimov et al, 33,34 we had compared the split-step method (SSM) with the iteration method propose by us, and we have shown that application of the SSM possesses some disadvantages for the absorption OB problem solution in 2D case. Among them, we stressed the absence of the conservatism property and rounding error accumulation. It means the asymptotic stability property violation.
Earlier, in Trofimov et al, 35 we have investigated the 2D problem of absorption OB, taking into account the longitudinal and transverse diffraction effects. In this case, the mathematical model of the problem contains the Schrödinger equation with respect to the envelope (complex amplitude, slowly varying in time only) of a wave packet instead of PDE concerning the laser pulse intensity. There, we demonstrated existence and uniqueness of the solution for the differential problem in 2D case, positiveness of the free-charge concentrations in differential and difference cases, uniform boundedness of the mesh functions at the iterations, the convergence of the developed two-stage iteration process, and existence of the difference problem solution.
So, in Trofimov et al, 33 we discussed a 2D model of a special type of the absorption OB-field OB, which corresponds to the semiconductor absorption coefficient nonlinear dependence on the electric field potential. Two rules of the electric field potential re-normalization were proposed to avoid the Neumann problem solution multiplicity. Let us note that at computer simulation, we also pay attention to the uniqueness of the Neumann problem solution.
It should be stressed that in our previous papers, we did not discuss the direct numerical methods for the solution of the Poisson equation with the Neumann BCs. Therefore, in the present paper we pay a special attention to a possibility of using this method for the solution of the Poisson equation involving in a set of the nonlinear equations, describing a process of semiconductor plasma generation, and discuss a correlation between conservatism property of the problem under consideration and solvability condition for the Neumann problem.
The paper is organized in the following way. The mathematical model of the problem is discussed in Section 2. In Section 3, we develop and investigate the conservative FDS for the problem solving and construct the original threestage iteration process (TSIP) for its resolvability. The computer simulation results are presented and discussed in Section 4.
It should be stressed that many other processes, such as the problem of cell chemotaxis, 36,37 for example, or time-depended advection-diffusion-reaction problems, are described by the similar PDEs. So the developed numerical approach could also be applied for solving problems belonging to various scientific areas and have a wide application.

| PROBLEM STATEMENT
We consider the problem in the parallelepiped domain The process of semiconductor plasma generation induced by a femtosecond laser pulse is described by the following set of 3D dimensionless nonlinear PDEs 1,2 : ∂I ∂z Variables x, y, and z are the dimensionless spatial coordinates, and parameters L x , L y , and L z denote their maximal values, respectively. Variable t denotes a dimensionless time. This set of equations consists of nonlinear nonstationary equation concerning a free electron concentration n(x,y,z,t) in the conduction band of a semiconductor, the equation concerning an ionized donor concentration N(x,y,z,t). The Poisson equation is written with respect to a function φ(x,y,z, t), which describes a dimensionless electric field potential. The equation, describing the laser radiation intensity I(x,y,z, t) evolution at the pulse propagation along the z-axis, is written without taking into account the diffraction effects. This is valid if the beam diffraction length is many times greater than the semiconductor length. The electron diffusion coefficients D x , D y , and D z and the electron mobility coefficients μ x , μ y , and μ z are nonnegative parameters. Parameter γ depends on the maximal concentration of free charged particles, in particular. Parameter δ 0 denotes a maximal value of the absorption coefficient for laser energy.
The BCs for Equations (1) to (4) are written below: ∂φ ∂x ∂n ∂x −μ x n ∂φ ∂x They correspond to the absence of electric current through the semiconductor faces and the external electric field absence.
Incident laser pulse initial distribution is stated in the form of Gaussian beam with time-dependent pulse shape: From (7), we see that the pulse intensity at the beam center (x = L x /2,y = L y /2 ) tends to unite with time increasing.
The initial conditions for the functions are written in the following manner: where parameter n eq denotes an equilibrium value of the free charged carrier concentration. The generation G(n,N,I) and the recombination R(n,N) of free-charged particles are described by the following functions: where parameter τ R characterizes a free electron recombination time and q 0 is a maximal intensity of the incident laser pulse.
It should be stressed that because of the incident intensity spatial distribution (7) and absence of the external electric field, the problem solution is a symmetrical one concerning the laser beam center along x-and y-coordinates. This fact is an important feature for estimating the FDS's accuracy.
A nonlinear dependence of the absorption coefficient δ(n,N) from the semiconductor characteristics can be approximated in different ways depending on physical experiment conditions. In this paper, we consider its approximation, which corresponds to the concentration absorption OB: The charge conservation law is valid for the problem under consideration if the charge flow through the semiconductor boundaries is absent and we suppose that a semiconductor is electrically neutral before the laser pulse action: As it is well-known, [17][18][19][20] the Neumann problem is solvable if the following equality satisfied: where dσ denotes the surface element. In the case under consideration, the condition (12) is written in the following form: where v is one of the spatial coordinates x, y, z and Therefore, at the statement of the uniform BCs (5), the condition (13) coincides with the charge conservation law (11). Thus, a property of the numerical method conservatism is a key point for the problem solution.

| Conservative FDS
To solve numerically the initial-boundary problems (1)(2)(3)(4)(5)(6)(7)(8), the uniform mesh in time and space coordinates can be defined as follows: where h x ,h y ,h z , and τ are mesh steps on spatial coordinates and on time respectively and P x ,P y ,P z , and P t are number of meshes nodes. Let us define the mesh functions n h ,N h , and φ h on the meshΩ as follows: A mesh function I h is defined on the meshΩ 00 shifted on spatial coordinate z and on time coordinate that is caused by the requirement of the FDS second-order approximation achievement for the equation describing the intensity evolution. For brevity, we use below the following index-free notations: where f is one of the mesh functions n h ,N h , and φ h , and the following ones concerning I h : At FDS developing, we also introduce some useful notations: The first and the second difference derivatives are defined in a standard way, and they are notated as follows: 3D Laplace difference operator is stated also in the standard way: For difference equations notation brevity, we introduce the finite-difference operators: and the similar operators on y-coordinate and z-coordinate.
Using the Crank-Nicolson scheme, we write down the symmetric FDS for the problem (1-8): The BCs are approximated as follows: Intensity initial distribution is approximated in the following way: Initial distribution for other functions is stated similarly as for differential problem: Theorem 1 is valid for the developed FDS: The FDS (22)(23)(24)(25) approximated the initial-boundary problem (1)(2)(3)(4) with the second order on spatial coordinates and on time coordinate in inner mesh nodes concerning the point (x i , y j , z k , t m +0.5τ) on sufficient smooth solution of the differential problem.
Let us note that using the standard differential derivatives expansion in a Taylor series on enough smooth solution of the problem (1) -(4), it is easy to prove this theorem. For brevity, we omit the proof of Theorem 1. Implementation of the difference analog for the conservation law (11) is very important for the problem numerical solution, and it means a conservatism property of the FDS. So we follow the difference invariant at the FDS construction. For its numerical representation, we use the trapezoid rule that possesses the second order of approximation: It should be stressed that we compute the invariant (29) in inner mesh nodes because of an optical radiation presence at two boundaries of the domain on z-coordinate. If the conservatism property does not violate due to the computational error accumulation at computation in a long-time interval, then the FDS possesses the asymptotic stability property. Let us note that the correlation between the FDS conservatism property (29) and the Neumann problem solvability condition (13) will be discussed below in Section 3.5.

| Proof of the FDS (22-27) conservatism.
It should be stressed that the BCs (26) approximate the BCs (5-6) with the first order of accuracy. This is caused by the requirement of the

Theorem 2 The FDS (22-27) is a conservative one at the BCs approximation (26).
Proof. Subtracting (23) from (22), we obtain the following expression: Therefore, Q(t m ) evolution is governed by the following equation: Let us consider the first term entering the sum above. We see that the following equalities Thus, the first term of Equation (30) right part can be transformed as follows: Substituting the BCs (26) in (31), we can state that for any mesh node t m > 0. In similar way, it is easy to show that other terms of (30) are also equal to 0. Taking into account the initial condition (28), it is obviously that Q(t 0 ) = 0. So we obtain that Thus, the charge conservation law (29) is valid, and the FDS (22-27) is a conservative one.

| Three-stage iteration method for FDS realization
There are various numerical methods for solving the problem similar to Equations (1) to (8). As it was mentioned above, the splitting methods are widely used. But these methods possess some disadvantages: for example, the conservatism property does not realize. So to avoid this, we develop the original TSIP for the conservative FDS (22)(23)(24)(25)(26)(27) resolvability. Below, we apply both methods for the problem solution and compare obtained computer simulation results.
The key difference between the SSM and the iteration method is demonstrated in Figure 1, where the function f denotes one of the mesh functions, which were introduced above: n h ,N h ,φ h ,I h .f and f and f are the functions on the upper time layer and on additional time layers. At using TSIP, we provide a computation of the problem solution on the upper time layer using function values obtained on the previous time layer and the previous iterations. In the case of using the SSM for the problem solution, the additional time layers are introduced, and we have to solve the difference equations on these additional layers before we obtain the solution on the upper time layer.
We notice that the upper indexes s, s + 1, and s + 2 over the mesh functions mean the iteration number. Let us write the stages of TSIP separately with corresponding BCs. Below, we use the following notation for operator Z 0:5 Þ and the similar notations for operators on other spatial coordinates. The first stage of the iteration process is as follows: At this iteration stage, the BCs (38) are used at solving the Eequation (35) by the meaning of the Thomas algorithm. The BCs (39) and (40) are used for computing the values of the functions on the domain boundary planes (0 ≤ x ≤ L x ; 0 ≤ z ≤ L z ) and (0 ≤ x ≤ L x ; 0 ≤ y ≤ L y ), respectively. The similar algorithm of the BCs using is applying for the second and the third iteration stages.
The second stage is as follows: The third stage is as follows: , For brevity, we will name all three stages (34-51) of the iteration process as its full stage. The iteration process is stopped if the corresponding convergence estimation (criterion of the relative error) is fulfilled:n As follows from (52), we compare the solutions for the full stage of the iteration process. If one of these inequalities are not valid for any function at any mesh node, then we repeat computation. As initial approach for the iteration process (34-51), we use the mesh functions values, obtained on the previous time layer:

Theorem 3 FDS (22-27) is a conservative one on each of the iteration stages if the BCs are approximated with the first order by (38) to (40), (45), and (50).
Proof. The technique of Theorem 3 proof is similar to Theorem 2. To proof the occurrence of the conservatism property on each of the iterations, we should consider the following differences: As it was shown above when Theorem 2 was proved, the invariant Q(t m ) on the lower time layer is equal to 0. So we should estimate the invariantsQ Þon the upper time layer: At stating of the zero-value BCs and their approximation by (38) to (40), (45), and (50), the right parts of equalities (55) to (57) are equal to zero. Thus, the conservatism property is valid for each stage of TSIP.
Notice, if a semiconductor is placed into the external electric field, which results in the nonuniform BCs for the problem under consideration, the conservation property validation on each of the TSIP stages is a more complicated problem, and it requires an additional investigation. It should also be stressed that the FDS conservatism on each of the iterations is a critical property for the asymptotic stability of the FDS.

| SSM for the problem solution
As it is well-known, the SSMs are widely used for a solution of multidimensional nonlinear problems. So to estimate the efficiency of TSIP, we compare the computer simulation results obtained by both methods. For this purpose, we write down a FDS based on the SSM for the problem (1-8) by introducing the additional time layers: (m+1/3)τ and (m +2/3)τ (see Figure 1). So to carry out a computation on the upper time layer (m+1)τ, it is necessary to solve the set of equations on these additional layers. Thus, we introduce the following time meshes: We also introduce additional notations: Below, the SSM for the solving of the problem (1-8) at additional time layers (m+1/3)τ and (m+2/3)τ and at upper time layer (m+1)τ is written in the form: Λφ for computation of the problem solution at the first additional time layer; Λφ for computation of the problem solution at the second additional time layer;

| Numerical methods for Poisson equation solution
At the system (1-8) investigation, we have to solve the 3D Poisson equation concerning the electric field potential with Neumann BCs. As the process under consideration is a nonstationary one, and this process is studied in a rather big time interval, then both in the case of TSIP using and in the case of SSM using, we need to find the numerical solution of Equation (1) at each of the time layers and on each of the iteration stages or on each of the additional time layers. In this regard, we should pay special attention to the computational efficiency of a method for the Poisson equation solution. As was noted above, the right part of the Poisson equation depends on the functions, which are computed with some computational error, and this error could accumulate and increase during computations. This can lead to the Neumann problem solvability condition violation. It is also should be noted that the accuracy of the electric field potential calculation is extremely important because of the nonlinear feedback presence between the charged particle concentrations and the electric field potential.
In our investigation, we consider two approaches for the Poisson equation solution with Neumann BCs: the direct method based on the FFT 21-23 and the iteration method on the base of a split-step scheme. [25][26][27][28] In the case of the direct method application, we have to solve the difference Poisson equation (22) for which the difference analog of the solvability condition (13) should be valid. Taking into account the uniform BCs (26), this condition coincides with the difference conservation law (29) and takes the form: Obviously, the conditions (76) could be carried out only with some computational accuracy because of rounding error accumulation despite the FDS conservatism property. Thus, at solving the nonstationary problem under consideration, we meet the solvability condition violation at using the direct method for the Poisson equation solution, which could lead to the distortion of the problem solution.
At using the iteration process, based on the stabilization method, the difference equation concerning electric field potential is written in the form: which is similar to the difference approximation of the Helmholtz equation It is well-known that the Neumann problem for the elliptic equation (78) is uniquely solvable if η > 0. Therefore, Equation (77) has a unique solution as τ is a positive constant. That is why we use the iteration method and demonstrate its advantages in comparison with the direct method using. Thus, we use the additional iteration process for the Poisson equation solution: We apply the Thomas algorithm to solve each of the difference equations (79-81). The convergence assessment for the iteration process (79-81) is given by the inequality which is the estimation for the difference Poisson equation residual (we discussed a choice of the convergence criterion in Trofimov et al 34 ). If the criterion (82) is satisfied for all mesh nodes, then the additional iteration process is stopped, and the function F b+1 represents the required electric field potential. As the solution of a problem with the Neumann BCs could be found only up to a constant, it can lead to the multiplicity of the problem (1-8) solutions. To avoid this uncertainty and to achieve the problem solution uniqueness, we provide re-normalization of the electric field potential values by means of the following formula: The iteration method (79-81) has the second order of accuracy concerning the iteration parameter τ.

| COMPUTER SIMULATIONS RESULTS
In this section, we discuss the computer simulation results, provided with different mesh steps to compare the efficiency of investigated numerical methods. We estimate the properties of the methods by using rather big mesh steps and accuracy of the conservation law preservation. First, we consider the mathematical model corresponding to zero-value electron mobility. In this case, nonlinear feedback between the functions of the system is absent. Therefore, the concentration of free electron does not depend on the electric field strength. In Figure 2 (3) concerning the free electron concentration. Therefore, the increasing of computational errors does not occur. However, the method for the problem (1-8) solution plays an essential role. As it is seen from Figure 2, the accuracy of invariant conservation is six orders higher at using TSIP than in the case of using SSM. The maximum Q(t m ) value at using SSM ( Figure 2, red dotted line) is about 2.3 Á 10 −6 , and it is reached at the time point corresponding to the time moment of the explosive absorption appearance in the OB system and formation of a contrast spatial-time structure.

| The case of electron mobility presence
Next computations are provided with nonzero electron mobility. In this case, a method choice for the Poisson equation solution significantly influences the accuracy of the invariant conservation. As can be seen in Figure 3, at using the iterative process (79-81), we obtain the same invariant evolution as the one presented in Figure 2. We want to emphasize the extremely high accuracy of the invariant conservation at using TSIP for the problem (1-8) solving, which confirms the conservativeness of this method. Moreover, at using FFT algorithm for the Poisson equation solution, another situation is observed (Figure 4). The invariant Q(t m ) grows with a dependence close to linear. This can lead to an essential distortion of the numerical solution at providing computations in a long-time interval. Decreasing the mesh time step τ by order does not result in significant improvement of the invariant conservation accuracy ( Figure 4B). It should be stressed that at solving the system (1-8) by TSIP, the invariant Q(t m ) conservation accuracy is better on seven orders if we solve the Poisson equation by using the iteration process (79-81) than if we use FFT algorithm. This occurs because of the accumulation of computational errors in the right part of the Poisson equation and (as a consequence) the violation of the conservatism property. Therefore, the Neumann problem solvability condition is also violated. Thus, it can be concluded that for the solution of the Neumann problem for the Poisson equation, involving in a set of time-depended nonlinear PDEs, the iteration method application is much preferable than using of the direct method, because the last one results in the conservatism property violation.

| Comparison of methods efficiency for the problem (1-8) on the base of functions norms estimation
Below, we compare the efficiency of using TSIP and SSM for the problem (1-8) solution. At computer simulation, we use the iteration process (79-81) for the Poisson equation (1) solving. We choose the electron mobility being equal to unity: μ x = μ y = μ z = 1. We estimate the changes in the problem solution at the mesh step τ decreasing. The norms of the difference between the free electron concentrations computed with different τ are presented in Figures 5-7. As is seen in Figure 5, at using TSIP for the problem solution, the mesh time step decreasing by order of magnitude does not lead to significant solution changing: kn(τ 1 ) − n(τ 2 )k C has the order of about 10 −5 , and n τ 1 ð Þ−n τ 2 ð Þ k k L 2 has the order of about 10 −7 , where τ 1 = 10 −3 and τ 2 = 10 −4 . Let us note that the norms are computed by formulas F I G U R E 5 Difference of norms for the free electron concentration computed using threestage iteration process (TSIP) with τ 1 = 10 −3 and τ 2 = 10 −4 F I G U R E 6 Difference of norms for the free electron concentration computed using split-step method (SSM) with τ 1 = 10 −3 and τ 2 = 10 −4 If we use SSM for the computations, then these norms have the orders of about 10 −2 and 10 −4 , respectively ( Figure 6).
Let us analyze the computer simulation results, presented in Figure 7. It demonstrates the norm of the difference between solutions, computed by using TSIP with τ 1 = 10 −3 and by using SSM with different mesh steps: τ 2 = 10 −3 or τ 2 = 10 −4 . As one can see, the solutions, obtained on the base of TSIP and SSM, differ strongly if we use the equal time steps τ 1 = τ 2 = 10 −3 (Figure 7, black solid line). If we decrease the time step up to τ 2 = 10 −4 at using SSM, the solution converges to the solution obtained by TSIP at using τ 1 = 10 −3 (Figure 7, red dashed line). As it is also seen, the norm values increase with time. This is caused by the accumulation of the computational error and a violation of the conservativeness property in the case of using SSM. Thus, TSIP is an effective method for a solution of the problem under consideration. It is a conservative one and allows us to provide computations with a rather big mesh step on time coordinate without losing the solution accuracy, which is a principal advantage for the modeling of multidimensional problems.

| Formation of contrast spatial-temporal structures in semiconductor
Finally, in this paragraph, the formation of complicated contrast spatial-temporal structures is demonstrated in Figure 8. The computations are provided by using TSIP. The Poisson equation (1) with corresponding BCs is solved by using the iteration process (79-81). To investigate such periodic (quasi-stationary) modes, it is necessary to carry out computer simulations in a long-time interval (t ≥ 1000). Therefore, the requirement of the numerical method conservatism and its asymptotic stability is a crucial point for computer modeling of the problem under consideration. As it was shown above, TSIP is an effective method for the absorption OB problem numerical solution. The asymptotic stability of TSIP is confirmed by maintaining the symmetry of the solution, which corresponds to the physical meaning of the problem in the case of uniform BCs. The high accuracy of invariant Q(t m ) preservation occurs even at providing the computations in a long-time interval. F I G U R E 7 Difference of norms for the free electron concentration computed using threestage iteration process (TSIP) with τ 1 = 10 −3 and split-step method (SSM) with τ 2 = 10 −3 (black solid line) and using TSIP with τ 1 = 10 −3 and SSM with τ 2 = 10 −4 (red dotted line) [Colour figure can be viewed at wileyonlinelibrary.com]

| CONCLUSIONS
In the present paper, the numerical methods for solving the 3D problem of laser pulse interaction with a semiconductor are discussed. This interaction is described by a set of PDEs involving linear elliptic equation with zero-value Neumann BCs, the reaction-diffusion equation, containing a term like a thermo-diffusion term, and two other nonlinear equations.
We showed that despite a wide application of the direct method (FFT) for solving the Poisson equation with zero-value Neumann BCs, this method could not be used in the case under consideration because of the solvability condition violation for this equation. It arises from the charge conservation law violation because of the presence of both round-off errors at the computer simulation and the feedback between changing of the functions involved in the set of the nonlinear equations. Therefore, a numerical method for the Poisson equation solution with zero-value Neumann BCs must be based on the iteration method (SSM) application. In this case, the solvability condition is satisfied. For computer simulation of the 3D problem of femtosecond laser pulse interaction with a semiconductor, the conservative FDS is constructed, and the TSIP for its resolvability is proposed. This iteration process allows us to develop the numerical method economical like an SSM. Moreover, the proposed FDS possesses an asymptotic stability property. The FDS conservatism property is proved, and it is also confirmed by the computer simulation results. Validation of this property is crucial for the solution of the problem under consideration because it coincides with the condition of Neumann problem solvability.
In our opinion, the developed TSIP is an effective tool for the solution of nonlinear PDEs. The numerical approach developed in the paper can find applications in other fields of science.