Abstract

In this study, we introduce an innovative approach to solving stochastic equations in two and three dimensions, leveraging a time-splitting strategy. Our method combines radial basis function (RBF) spatial discretization with the Crank–Nicolson scheme and the local one-dimensional (LOD) method for temporal approximation. To navigate the probabilistic space inherent in these equations, we employ the Monte Carlo method, providing accurate estimates for expectations and variations. We apply our approach to tackle challenging problems, including two-dimensional convection-diffusion and Burgers’ equations, resulting in reduced computational and memory requirements. Through rigorous testing against diverse problem sets, our methodology demonstrates efficiency and reliability, underscoring its potential as a valuable tool in solving complex multidimensional stochastic equations. We have validated the method’s stability and showcased its convergence as the number of collocation points increases. These findings serve as compelling evidence of the suggested method’s convergence properties.

1. Introduction

The realm of mathematical modeling is often confronted with the inherent uncertainty introduced by perturbations in various physical phenomena. To establish dependable models capable of capturing this uncertainty, the utilization of stochastic partial differential equations (SPDEs) becomes imperative [1–3]. The emergence of SPDEs as a captivating research domain in recent years has garnered significant attention from scholars worldwide.

While the literature on SPDEs is extensive, the scarcity of exact solutions necessitates the use of numerical techniques for their resolution [4–8]. Researchers such as Yoo [9], Yan [10], Nouy [11], and Ye [12] have applied diverse numerical methods, including finite difference, finite element, and kernel-based collocation methods, to address the challenges posed by stochastic PDEs.

Among the array of numerical techniques available, mesh-free methods, which have gained substantial popularity, hold promise for tackling these complex equations [13]. Specifically, the collocation method based on radial basis functions (RBFs), belonging to the class of mesh-free methods, has been widely adopted by researchers for both deterministic and stochastic PDEs [14–18]. In this study, our main focus is on solving stochastic diffusion equations using the collocation method based on radial basis functions (RBFs). The RBF method has a distinct advantage as it does not rely on mesh discretization and can be applied to unstructured node sets. This makes it particularly suitable for solving problems in complex and irregular domains, where geometric details may be intricate. Table 1 provides definitions for the most widely used RBFs, where the Euclidean norm is denoted by , and the shape parameter controls the influence of the RBF.

The computational cost associated with solving two- or three-dimensional partial differential equations (PDEs) is known to be significant. However, this challenge can be effectively addressed through the utilization of time-splitting methods. These methods aim to reduce the dimensionality of the algebraic system involved in the numerical scheme. By decomposing complex time-dependent problems into simpler subsets, the numerical methods can be solved independently and efficiently. Then, the existing numerical methods can be solved separately and shortly. Notably, recent studies have explored the application of time-splitting approaches for multidimensional equations [19, 20]. In this particular investigation, we employ the local one-dimensional (LOD) time-splitting technique to tackle the two- and three-dimensional versions of stochastic diffusion equations [21]. In the subsequent subsection, a concise overview of stochastic PDEs will be presented.

1.1. Stochastic Diffusion Equations

Consider the stochastic advection-diffusion equation as a popular stochastic problem. Let be a probability space, where is the space of basic outcomes, is a -algebra related to , and is a measure defined on . Many phenomena in biology, chemistry, physics, and population dynamics can be expressed in the format of a stochastic diffusion equation with the following form:where is a bounded domain and is a constant. is a two-dimensional random noise related to the Brownian motion , which is assumed that it depends on time and space with zero mean and spatial covariance function given by

Equation (1) represents a well-known problem in industry and mathematics that is a focal point for numerous researchers who have endeavored to find its optimal solutions or have utilized this set of equations as a means to showcase the effectiveness of their proposed methodologies. Discussions have already unfolded concerning the existence and uniqueness of solutions within this class of equations [22, 23]. To describe the white noise different forms can be used of the covariance function. In this technique, we define it in a way that is proportional to the chosen radial basis function.

The structure of this paper is outlined as follows: Section 2 provides a concise overview of the implementation of the radial basis function (RBF) method, while Section 3 introduces the local one-dimensional (LOD) method for addressing time-dependent problems. In Section 4, we present the general form of stochastic diffusion equations. To assess the efficacy of the LOD-RBF method, we employ illustrative examples and analyze its performance. Section 5 concludes this paper by summarizing the key findings and highlighting the main conclusions derived from our study.

2. Collocation Method Based on Radial Basis Functions

In this section, we elucidate the process of approximating numerical solutions using the radial basis function (RBF) method for one-dimensional problems. Subsequently, this procedure can be extended to encompass two-and three-dimensional problems.

To illustrate, let us consider a stochastic diffusion equation expressed in the following form:where is the diffusion coefficient and is a nonlinear function. Furthermore, are given functions, and is the computational domain.

For time discretization, we apply the Crank-Nicolson scheme between two consecutive time steps. Let be the time step size

For , applying a forward finite difference scheme for (3) leads to the following equation:where

Here, the random variable , has the following mean and covariance function:

After time discretization, the Wiener process transforms into a Gaussian field. Then, in order to present a space discretization for (6), we apply a collocation method based on the radial basis function at fixed collocation points. We choose the collocation pointswhere for belongs to the domain and for belongs to the boundary. Considering as the Euclidean norm, the RBF method strives to approximate solution using a linear combination of RBFs at the the collocation points. Thus, we define the following estimate:where is a radial basis function that is defined by . The coefficients , are unknown, and they should be determined at each time step.

The approximation (10) can be written at points as follows:which can be written in the following matrix form:where

Considering the above notations, for , the RBF collocation method for the discretization of (6) reads as: Find such thatwhere

Moreover, for the boundary conditions we have

The equations result in the following nonlinear system:where is defined as follows:and the matrices and and the right-hand side vector can be obtained from (14)–(18) as follows:

In order to solve the nonlinear system (19), we use a Newton’s method.

3. The Locally One-Dimensional Time-Splitting Technique

In this study, we employ the local one-dimensional approach. Specifically, we examine the diffusion equation presented as follows:

In the domain and certain boundary and initial conditions. are constant fluid velocities and are the dispersion coefficients in the x-direction and y-direction, respectively.

We split the above equation into two one-dimensional equations as follows:

Each of the above equations should be solved over half of the time step to obtain the approximated solution for the 2D advection-diffusion problem. These equations are easily solved using the schemes developed for the 1D ones. In addition, instead of solving (23) and (24) in each half-time step, the following equations can be solved over a full-time step:

To solve (25) and (26), we can apply any methods used for solving the one-dimensional advection-diffusion equation. In the same way, by applying the LOD approach to (3), we have two one-dimension equations as follows:where

In order to solve the nonlinear system (27) and(28), we use a Newton’s method. We wish to emphasize that our choice for the initial value of the Newton method was derived from a coefficient within the problem’s initial conditions.

Before delving into the subsequent section, where we will unveil our numerical findings, we embark on an analysis of the method’s convergence. To accomplish this, we turn your attention to the following theorem.

Theorem 1. The proposed approach RBF-LOD is stable.

Proof. While maintaining the problem’s generality, we focus on the homogeneous case, noting that the same procedures can be readily adapted for the heterogeneous case with only minor adjustments. Hence, we proceed under the assumption thatAlso considerwhere represents the angular part of the points . We can rewrite in the following form:whereand is a radial basis function. We haveNow, upon substituting (34) into the modified homogeneous (3), we obtainLet us rephrase (35) for each .The resulting system can be represented in matrix form as follows:whereUnder the assumption of the following, we will obtain:Equation (37) can be expressed as follows:Since A is a positive definite matrix, it is invertible, enabling us to express:Given that the term represents a coefficient matrix with a condition number in each iteration that does not exceed a certain value, denoted as , we can state the following:Now, employing an iterative numerical method, denoted as , we can decompose the equation above into discrete time steps as follows:where denotes the number of iterations in the numerical method. In (43), the matrix is influenced by the matrix . The relationship described in (43) is Lax-stable ifIf is a normal matrix, establishing the stability of the method requires demonstrating that the eigenvalues of fall within the stability domain of the method .

Lemma 2. A, B, and C are commutative matrices.

Proof. For proof refer to [24]

Lemma 3. and are skew-adjoint.

Proof. From (38), A and B are skew-adjoint. Furthermore, since A and B commute by Lemma 2, so too do and B. Hence,where denotes the complex conjugate transpose. Using the same methodology, the proposition is also proven for the matrix product .

Lemma 4. and are Hermitian matrices, and consequently, they are Normal matrices.

Proof. For proof refer to [24]

Theorem 5. The relation (43) is Lax-stable.

Proof. Matrix is skew-adjoint so it is a Hermitian matrix. As stipulated in Lemma 4, due to the Hermitian nature of matrix , it also qualifies as a normal matrix, ensuring that its eigenvalues are real. Therefore, for equation (43) to be Lax-stable, it suffices for the eigenvalues of matrix to reside within the stability domain of method . By making a judicious choice of method , such as the Crank–Nicolson method, which boasts unconditional stability with a domain encompassing all real numbers, we can confidently assert that equation (43) indeed is Lax-stable.
In light of the preceding explanation, we can confidently assert that the proposed method is stable.

4. Numerical Results

In this section, we showcase the numerical results obtained through the utilization of the LOD-RBF method. Our primary objective is to obtain accurate numerical solutions for stochastic diffusion problems. We demonstrate the effectiveness of the proposed method by solving various stochastic diffusion equations in both two- and three-dimensional space. We employ a Monte Carlo simulation for the solution samples for , and estimate the expectation of the solution at the final time bythat is called the mean value. Also, we compute the variance and the standard deviation for all test problems presented that are defined as follows:knowing them as are two next of the statistical moments is momentous. To assess the accuracy of the obtained numerical solutions, we utilize two commonly employed error measures: the root-mean-square error (RMSE) and the maximum error . These metrics provide valuable insights into the quality and precision of the numerical results. To define these errors, first, we introduce as the expected error function at time in the following formulate:where and are the exact and expectation of numerical solutions at the collocation point and time , respectively. The discrete errors are defined as follows:

Example 1. Convection-diffusion models play a fundamental role in describing the transport phenomena encountered in various physical, chemical, and biological processes [25]. In this study, we focus on the two-dimensional stochastic convection-diffusion equation (2D) with Dirichlet boundary conditions. This equation serves as a key framework for analyzing the intricate interplay between convection and diffusion processes in diverse fields of study. Consider the following stochastic convection-diffusion equation:where is the Reynolds number that depends on the specific problem and the fluid flow conditions involved. The exact solution of the deterministic type of problem isHere, we approximate the numerical solutions by using inverse multiquadric (IMQ) functions. Concerning the covariance function, due to the IMQ basis function, we consider it in the following form:Also, our selected optimal shape parameter is which we determined by trial and error. Other required assumptions and constants are as follows:We get all the results by using realizations of the Monte Carlo method. In Figure 1, the exact solution, the mean solution, and the error are shown. Furthermore, the standard deviation and variance are shown in Figure 2.
The CPU times for implementing this approach in each example are provided. The computations were performed on a computer equipped with a 16-core neural engine and 16 GB of memory using MATLAB 2021b. According to the aforementioned constants, the CPU time for implementing this approach in the specified example is 0.71 seconds. In Table 2, the efficacy of the number of collocation points M on the errors (RMSE and ) is investigated. As expected, the increase of collocation points leads to the decrease of error.
It is evident that employing a smaller time step size leads to improved accuracy in the obtained results. This observation is substantiated by the data presented in Table 3, which clearly demonstrates a reduction in errors as the time step size decreases.
In order to assess the stability of the method, we conducted a study, considering the suggested optimal shape parameter as the input data. Our aim was to examine whether the condition number of the matrix, resulting from solving the system of linear equations, remains relatively unchanged in near the selected optimal shape parameter. The results supporting this proposition are presented in Table 4. Based on the findings obtained from these experiments, we conclude that employing ensures both stability and high accuracy in the numerical solutions.
Furthermore, it can be intuitively demonstrated that the proposed method for this class of equations is convergent. To substantiate this claim, we illustrate how the error decreases as the number of collocation points M increases.
As vividly illustrated in Figure 3, it is evident that the error decreases with the increasing number of collocation points. Consequently, this observation serves as compelling evidence of the method’s convergence. In this particular example, we have achieved a convergence of orders 3 and 4.

Figure 1 

Exact solution, mean solution and, error for , and .