On a Method for Solving of Multidimensional Equations of Mathematical Physics

. We consider linear multidimensional evolutionary equations or the linear part of nonlinear ones. The complex structure and the presence of terms with di ﬀ erent physical sense requires the coordinate splitting to be preceded by splitting by physical factors (processes). In contrast to the coordinate splitting this kind of splitting can be exact in some nodes and in the intervals it can be controlled. The method in question is developed in 80s of 20th century by G. Marchuk [1] and till now it is applied very successfully for solving various problems in ecology, air and water pollution, di ﬀ usion, etc. Here we show that this method is relevant for studying of propagating ultrashort localized pulses in nonlinear waveguides as well.


Introduction
In this paper we aim to demonstrate that the splitting by physical factors is applicable and can be efficient for solving multidimensional evolutionary problems in the optics. These equations possess very rich phenomenology, they contain nonlinear terms and do not admit analytic solutions being nonintegrable in the general case. Without loosing of generality we pay attention to an initial-value problem of (3+1)d equations of Schrödinger kind. We split the linear part of the differential operator and the initial conditions to two consequent Cauchy problems and by using a spectral analysis we proof that the original and the resulting splitting problems are equivalent in the ends of a given interval. For constant coefficients the above assertion holds in the whole interval. In opposite case for small enough step the approximation error varies in an acceptable range. The nonlinear terms are linearized by using of so-called inner iteration.
The method is applied successfully for numerical solving (3+1)d Schrödinger equation with sign-variable group velocity as well for (3+1)d amplitude equations (of envelope), when the propagation regimes of the light pulses may be ultra-short (femtosecond). The obtained results are reliable and give good predictions for the material quantities and dynamics of the light pulses.

(3+1)d Schrödinger Equation
Let us consider the Cauchy problem u(x, y, τ, z = 0) = g(x, y, τ), u| Γ = 0 where Γ is the boundary of 3d spatial region, g is a given function, τ = t − z v g is the local time in a moving frame with group velocity of the pulse v g , k is the wave number, β is a e-mail: mtod@tu-sofia.bg the dispersion, and z is the evolutionary coordinate of the pulse. Eq. (1) describes two different physical process: the first one is the diffraction of the pulse represented by the equation The second process relates to the dispersion and is governed by the equation When the initial function and the sought solution approach fast zero uniformly in all arguments there exist Fourier integrals whereĝ andû are the Fourier images of the functions g and u.
On the other hand the triple Fourier transform of Eq. (1) with respect to x, y, and τ gives a Cauchy problem of the ordinary differential equation which admits the exact solution Then the solution of the original problem, Eq. (1) has the following integral form

Fourier Analysis, Splitting by Physical Processes
Let us introduce an interval 0 ≤ z ≤ ξ and consider the following initial-value problems If function g decades on infinity then the Fourier-images of above problems exist. They admit exact solutionŝ After obvious algebraic manipulations we obtain For z = ξ the solutions (5) and (9) are equivalent for any ξ while for 0 ≤ z ≤ ξ the solution of the split problem is approximate. The latter means that in a discrete set of nodes z j , j = 0, 1, 2, ... one has to solve consequently the initial-value problems This is the gist of the method of splitting by physical processes applied for multidimensional problems.
By using Fourier analysis we showed that the original problem (1) is equivalent to solve consequently initialvalue problems (10) and (11) at a given set of nodes z j , j = 0, 1, 2, .... In the concrete case (1) is (2+1)d and needs one more splitting -this time coordinate and after that it can be solved numerically by alternative directions implicit (ADI) method, e.g. by the Peaceman-Rachford scheme or by the scheme of stabilizing correction [2]. In opposite, Eq. (11) is (1+1)d and can be solved direct by Crank-Nicolson scheme in complex arithmetic. The resulting three equations possess sparse matrices and can be easily treated by Thomas algorithm and its complex-valued modifications [3,4].
Splitting by physical processes is applied for more complicate Eq. (1) with cubic nonlinearity and an initial condition -function g(x, y, τ) ≡ exp . To get an implicit scheme about the nonlinear part also we use so-called internal iteration, which is convergent for small enough steps in z-direction [5,6].

(3+1)d Envelope Equation
Let consider one more Schrödinger-like equation -the envelope equation, which is (3+1)d dispersion equation with cubic, quintic and rational nonlinearities. The linear part of the equation is where ∆ = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 , a 1 , a 2 , a 3 are dimensionless real constants. For simplicity we will not comment their physical sense in this paper. The splitting by physical processes leads to following two initial-value problems The Fourier analysis results in consecutive solving in discrete set of nodes z j , j = 0, 1, 2, ... of initial-value problems u(x, y, τ, z = z j ) = u j+1 1 .
Eq. (15) requires a coordinate splitting, while Eq. (16) can be solved directly by complex-valued Crank-Nicolson scheme. Due to the higher dispersion the band matrix is four-diagonal. In order to use again Thomas algorithm modifications the mixed derivative in both problems should be put on the old z-stage. In this way the scheme of stabilizing correction for Eq. (15) yields i Here (·) stands for the values of the grid function on the z- The first equation is implicit in x-direction, while the second one -in y-direction. The fully fledged (3+1)d nonlinear Envelope equation is quite more complicate compared to (3+1)d nonlinear Schrödinger equation The last three terms describe the ionization processes, the cubic and quintic nonlinearities -the effects of self-steepening, the spatial self-focusing (self-defocusing), b 1 ,...,b 7 are real constant coefficients, ρ is the electron concentration, which is a dynamical quantity. The latter requires to solve Eq. (18) together with the kinetic equation ∂ρ ∂τ = t 1 (u)ρ+t 2 , t 1 -functional coefficient, t 2 -real constant.
(19) The nonlinearities can be linearized again by an inner iteration convergent for small enough z-steps. Eq. (19) being 1d for given values of set function u is treat by Crank-Nicolson scheme. An important advantage of the linear parts of the above splitting schemes is their absolute stability and full approximation of the differential operators.
The above considered equations are femtosecond (1fs = 10 −15 sec). In order to emphasize the significance of these equations we stress that they govern and are a tool for study of the spatio-temporal dynamics of the molecular and atomic processes.

Discussion and Conclusion
We implement the above splitting in finite differences for (3+1)d nonlinear Schrödinger equation and (3+1)d nonlinear envelope equation considering real regimes of pulse propagation. In this paper we do not intend to describe the physical results because they are considered thoroughly in our earlier papers [7], [8] and [9]. We aim only to motivate the splitting by physical processes combined with the inner iteration as a reliable method for solving of multidimensional wave and amplitude equations governing ultrashort regimes when exact solutions are impossible and the experimental observations and measurements are difficult and expensive.