Numerical Method for Stochastic Nonlinear Schr¨odinger Equation Driven by Multivariate Gaussian Measure: Algorithms and Applications

In this paper, we present a novel Galerkin spectral method for numerically solving the stochastic nonlinear Schr¨odinger (NLS) equation driven by multivariate Gaussian random variables, including the nonlinear term. Our approach involves deriving the tensor product of triple random orthogonal basis and random functions, which enables us to efectively handle the stochasticity and nonlinear term in the equation. We apply this newly proposed method to solve both one-and two-dimensional stochastic NLS equations, providing detailed analysis and comparing the results with Monte Carlo simulation. In addition, the proposed method is applied to the stochastic Ginzburg–Landau equation. Our method exhibits excellent performance in both spatial and random spaces, achieving spectral accuracy in the numerical solutions.


Introduction
Nonlinear Schrödinger equation is applied in many felds of physics, such as nonlinear modulation of collisionless plasma, nonlinear wave model, nonlinear optics, and logarithmic Schrödinger equation [1][2][3][4][5].In practice, some parameters in the physical equation may be uncertain or the model itself may have random disturbance and these disturbances may be related to each other.At present, the stochastic nonlinear Schrödinger equation studied are mostly driven by the Wiener process, white noise, or colored noise [6].It has been studied theoretically and numerically and takes the form where i � �� � − 1 √ , λ ∈ R, and W is a infnite-dimensional Wiener process.
Te polarization orientation of laser plasma can be characterized by initial conditions that exhibit both randomness and interrelationships.We can use the model as follows [7]: where x is the spatial variable and A + and A − are the clockwise and counterclockwise circularly polarized components, respectively.Te initial condition of A + and A − may have the Gaussian input as Te random variables α and β, which are correlated with each other, are often assumed to follow a multivariate Gaussian distribution in practical scenarios.Consequently, investigating the stochastic nonlinear Schrödinger equation driven by multivariate Gaussian random variables has become a prominent research topic.When solving this equation, the presence of the nonlinear term poses a challenging problem, regardless of whether it is driven by noise or a Gaussian random variable.
In [2], the NLS equation ( 1) is driven by the highdimensional colored Wiener process.In the Hamiltonian energy space, by implementing a scalar transformation on the noise, the uniqueness of the solution of the unknown function and the conditions for strong and weak solutions can be deduced theoretically [8].In [3], this paper focuses on the study of the nonlinear Schrödinger equation with a point-like source term.Te calculation of the nonlinear terms is closely linked to the energy norm.
In [9], the stochastic NLS equation (1) driven by white noise is numerically solved by symplectic and multisymplectic methods.Te calculation of the nonlinear term in the symplectic method involves satisfying the global Lipschitz condition, while in the multisymplectic method, it is determined by ensuring the discrete charge conservation law.Te two discrete schemes have good properties in longterm calculation, and both meet the discrete charge conservation principle [8,10].
In [11], this paper presents a proof that the solution to the stochastic nonlinear Schrödinger equation can be effectively approximated by the solution of a coupled splitting scheme.Numerical calculations show that both the evolution of charge and the exponential moments of energy reach strong convergence rates.In [12], this work employs the splitting approach to approximate the NLS equation and aims at demonstrating the strong and weak order of the splitting scheme.Convergence analysis is often used in the calculation of nonlinear terms, or the stochastic value of wave function at a certain time is directly used [3].
In this paper, we consider a stochastic nonlinear Schrödinger equation driven by multivariate Gaussian random variables of the following form: where a ∈ R, ∆ represents the Laplace operator, V( x → , ξ) is some known random function, and x → is the spatial vector.ξ is a multivariate Gaussian random variable.When the higher order numerical method is required, the stochastic Galerkin (SG) spectral method is favorable [13].Te numerical calculation of the nonlinear term and how to construct orthogonal basis is a challenging problem.In previous studies, many techniques have been proposed including the probability measure transformation method [14][15][16], domination method [14,[17][18][19], measureconsistent polynomial chaos expansion method [20,21], and Gauss-Schmidt orthogonalization method [22][23][24][25][26][27][28].
where α q   � α!/q!(α − q)!, and the other two parentheses of (6) do the same factorial operation.Te deduced tensor product of a double orthogonal basis and random functions is applicable when ζ is a multivariate Gaussian random variable.
Expanding on prior research, this paper aims at deducing the tensor product form of a triple random orthogonal basis (7) and random functions (8).We proposed a Galerkin spectral method for efectively solving the stochastic nonlinear Schrödinger (NLS) equation, which is driven by multivariate Gaussian random variables.In addition, we conduct numerical computations of the nonlinear term by representing it as a random function: where u(•), v(•), and w(•) are all random functions.Our main techniques include the following: (i) We construct a mapping transformation that connects multivariate Gaussian random variables to independent Gaussian random variables.(ii) Initially, we assume that the unknown function in the stochastic NLS equation can be expressed as a generalized polynomial chaos expansion.We then utilize the mapping transformation to convert it into an equation driven by independent Gaussian random variables.In addition, based on the derivations provided in Appendices A and B, we express the nonlinear term in the form of a generalized polynomial chaos expansion.Finally, we implement the stochastic Galerkin spectral method within the Gaussian measure space to solve the stochastic NLS equation.(iii) By solving the equations derived in the process, we obtain deterministic diferential equations for the coefcients of the expansion.
In addition, we employ the Monte Carlo method to simulate the stochastic NLS equation in one and two dimensions with diferent sample sizes.We then compare the results obtained using the proposed method with those obtained through the Monte Carlo simulation.
Te organization of this paper is as follows.In Section 2, we discuss the stochastic Galerkin spectral method for solving the NLS equation driven by multivariate Gaussian measure.Section 3 presents known conclusions in the feld.In Section 4, we deduce the tensor product forms of triple orthogonal basis and triple random functions.Section 5 demonstrates the application of our proposed method to solve one-and two-dimensional stochastic NLS equations.We analyze the results and compare them with those obtained using the Monte Carlo method.In addition, we apply the proposed method to the stochastic Ginzburg-Landau equation.Section 6 provides a summary of the proposed method and the results obtained.
(ii) Transform (4) into the equivalent variational formulation, for all orthogonal basis functions v, fnding ψ such that where the mathematical expectation E is calculated according to the joint probability density function of the measure ξ, and x → ∈ D ⊂ R d .(iii) Te weight function of v should correspond to the joint probability density function of ξ for orthogonal basis functions v i and v j satisfying the following equation for any i, j   ∈ J: So all of the v � φ κ , κ ∈ J N,K  , and the originally stochastic problem is reduced to the following fnite-dimensional and deterministic PDE problem, for all φ κ , fnding ψ N,K such that (iv) We use the second-order time splitting spectral method for the above deterministic problem.
We can construct the mapping relationship ξ ⟶ T η [31], which can potentially transform a dependent Journal of Mathematics multivariate Gaussian random vector ξ into a set of linearly independent Gaussian random variables η � (η 1 , η 2 , . . ., η K ).By the maps T, ξ � T − 1 (η), ( 9)-( 11) and ( 14) can be rewritten as follows: Te proposed method exhibits remarkable spectral accuracy convergence, particularly in its ability to achieve high precision in terms of the mean square error.To illustrate this concept, let us examine a scenario where both the spatial variable and the random variable are limited to a onedimensional context.Let u(x, t, ξ) and h(x, t, ξ) be the exact solution and the approximate solution, respectively [30,32].
where Q is a constant independent of M (truncation number), t is time, and m > 0 is a real constant depending on the smoothness of u in terms of ξ.
Given a fnite nonnegative integer N and K, we defne a truncated multi-indicates set and some set of operations [29].
Ten, any two polynomials from the set We know that the dimension number of the truncated polynomial space is To simplify the calculation, some important conclusions of orthogonal basis are established as follows.For any α, β ∈ J, Journal of Mathematics where B 1 is similar to (6) for multi-indicates situation; for more details of (25), see [29].Denote random functions u � u( x → , t, η) and v � v( x → , t, η) based on [29] uv � where θ 1 � α + β, q ≤ α ∧ β, and

Tensor Product of Triple Random Orthogonal Basis and Random Functions
In this section, we deduced the tensor product of a triple random orthogonal basis and random functions by relying on established conclusions.
In combination with Appendix A, for any α, β, γ ∈ J N,K , we have a multi-indicates situation as follows: Te proof of ( 28) and ( 29) can be seen in Appendix A.

Journal of Mathematics
Trough our numerical calculations, we have derived the following conclusions pertaining to the triple random unknown function.Suppose u, v, and w have the following SG chaos expansion: In combination with Appendix B, we have the following multi-indicates situation: where θ 2 � α + β + γ, p ≤ α ∧ γ, q ≤ α ∧ β, r ≤ β ∧ γ, and the multi-indicates situation Te proof of ( 32) and ( 33) can be seen in Appendix B.

Numerical Experiments
In this section, we test the accuracy of the proposed method by solving stochastic one-and two-dimensional NLS equations, respectively.All the presented stochastic equations are driven by bivariate random Gaussian variables.In addition, we assume that the unknown function have a periodic boundary condition in the spatial direction and will use the Fourier spectral method to approximate the unknown function ψ [33,34].
To minimize computational costs, we can derive the conclusions by utilizing equations ( 26) and (32), which are obtained from the application of equations ( 15)- (17).
Ten, ( 18) is equivalent to the Galerkin system as follows: To solve the Galerkin system (35) numerically, in the spatial and time direction of the Galerkin system, we used the Fourier spectral method and the second-order time splitting method, respectively.We used the pseudospectral method to compute the Fourier coefcients of  ψ κ− β+q  V β+q and  ψ α+p+q  ψ β+q+r  ψ γ+r+p [29].Te expectation, second-order moment, and variance of ψ( x → , t) can be calculated as follows: To demonstrate the accuracy of this method, we defne the errors in terms of expectation and variance for the real and imaginary parts, respectively.error and Var [•]  error correspond to the deviations in the real and imaginary parts of the expectation and variance, respectively.To emphasize the accuracy and efciency of the proposed method, we employ the Monte Carlo method to simulate the stochastic NLS equations in one and two dimensions using varying sample sizes.A comparative analysis is then performed between the Monte Carlo method and the proposed method..Consider the following stochastic one-dimensional NLS equation:
Equation ( 38) has the following analytical solution:

is taken in our calculation later).
Based on the maps T between ξ and η, Ten, f(ξ) and V(x, ξ) can be rewritten as follows: In Example 1, Tables 1-4 present a summary of our results, illustrating the errors in the real and imaginary parts of the expectation and variance of the random solution at diferent time points.On the other hand, Tables 5-8 display the results obtained through Monte Carlo simulations for equation (38).
Figure 1(a) displays the real and imaginary components of the expectation of the approximation solution.Figure 1(b) depicts the real and imaginary parts of the error in the expectation (top) and variance (bottom) for equation (38) at T � 2.

Stochastic 2D NLS Equation.
Consider the following stochastic two-dimensional NLS equation with the same random function f(ξ) as one-dimensional case:
Equation ( 43) has the following analytical solution: In Example 2, we present the summary of our results in Tables 9-12, which demonstrate the errors in the real and imaginary parts of the expectation and variance of the random solution at diferent time points.In addition, Tables 13-16 display the results obtained through Monte Carlo simulations for equation (43).
Figure 2(a) illustrates the real and imaginary components of the expectation of the approximation solution.By analyzing the data presented in Tables 1-4 and 9-12, along with Figures 1 and 2, it becomes evident that the proposed method achieves spectral accuracy in the context of the stochastic nonlinear Schrödinger equation with multivariate Gaussian measure, even in long-time scenarios.In comparison, Tables 5-8 and 13-16 indicate that Monte Carlo simulation does not exhibit high accuracy.Even when the accuracy occasionally reaches spectral accuracy, the computational cost remains signifcant and the overall effciency is low.

Figure 2 (
b) presents the real and imaginary parts of the error in the expectation (top) and variance (bottom) for equation (43) at T � 2, respectively.

Figure 1 :
Figure 1: (a) Plots of the expectation of approximation solution for equation (38) at T � 2. Te red and blue dashed lines represent the real and imaginary parts of the approximate solution, respectively.(b) Plots of the real and imaginary parts of the errors in the expectation (top) and variance (bottom), respectively.Te time step is ∆t � 0.0001.

Figure 2 :
Figure 2: (a) Plots of the expectation of approximation solution for equation (43) at T � 2. (b) Plots of the real and imaginary parts of the errors in the expectation (top) and variance (bottom), respectively.Te time step is ∆t � 0.000001.

Figure 3 :
Figure 3: Te expectation of amplitude of approximation solution for equation (46) at diferent time points for diferent parameters.Te time step is ∆t � 0.001.(a-d) Corresponding to the cases where α � 0.01, 0.1, 1.2, and 3, respectively.

Table 1 :
Te real part of the expectation error for the solution equation (38) at diferent time points by the proposed method and with time step ∆t � 0.0001.

Table 2 :
Te imaginary part of expectation error for the solution equation (38) at diferent time points by the proposed method and with time step ∆t � 0.0001.

Table 3 :
Te real part of the variance error for the solution equation (38) at diferent time points by the proposed method and with time step ∆t � 0.0001.

Table 4 :
Te imaginary part of the variance error for the solution equation (38) at diferent time points by the proposed method and with time step ∆t � 0.0001.

Table 5 :
Te real part of the expectation error for the solution equation (38) for diferent sample sizes and diferent time points by Monte Carlo simulation and with time step ∆t � 0.01.
S represents the sample size.

Table 6 :
Te imaginary part of expectation error for the solution equation (38) for diferent sample sizes and diferent time points by Monte Carlo simulation and with time step ∆t � 0.01.
S represents the sample size.

Table 7 :
Te real part of the variance error for the solution equation (38) for diferent sample sizes and diferent time points by Monte Carlo simulation and with time step ∆t � 0.01.

Table 8 :
Te imaginary part of the variance error for the solution equation (38) for diferent sample sizes and diferent time points by Monte Carlo simulation and with time step ∆t � 0.01.

Table 9 :
Te real part of expectation error for the solution equation (43) at diferent time points by the proposed method and with time step ∆t � 0.000001.

Table 10 :
Te imaginary part of expectation error for the solution equation (43) at diferent time points by the proposed method and with time step ∆t � 0.000001.

Table 11 :
Te real part of the variance error for the solution equation (43) at diferent time points by the proposed method and with time step ∆t � 0.000001.

Table 12 :
Te imaginary part of the variance error for the solution equation (43) at diferent time points by the proposed method and with time step ∆t � 0.000001.

Table 13 :
Te real part of expectation error for the solution equation (43) for diferent sample sizes and diferent time points by Monte Carlo simulation and with time step ∆t � 0.01.
S represents the sample size.

Table 14 :
Te imaginary part of the expectation error for the solution equation (43) for diferent sample sizes and diferent time points by Monte Carlo simulation and with time step ∆t � 0.01.
S represents the sample size.

Table 15 :
Te real part of the variance error for the solution equation (43) for diferent sample sizes and diferent time points by Monte Carlo simulation and with time step ∆t � 0.01.
S represents the sample size.

Table 16 :
Te imaginary part of the variance error for the solution equation (43) for diferent sample sizes and diferent time points by Monte Carlo simulation and with time step ∆t � 0.01.
S represents the sample size.