Perspective Chapter: Numerical Solutions for Modelling Complex Systems with Stochastic Differential and Partial Differential Equations (SDEs/SPDEs)

Abstract

Understanding phenomena ranging from biological processes to financial markets involves uncertainty. Stochastic Differential Equations (SDEs) and Stochastic Partial Differential Equations (SPDEs) serve as robust mathematical frameworks for modelling such systems. Given the stochastic influences within these models, comprehending the dynamics of complex systems becomes pivotal for accurately predicting system behaviour. However, traditional numerical techniques frequently encounter challenges in effectively addressing the intricacies and stochastic properties inherent in these equations. This chapter explores several numerical methods that offer streamlined and dependable solutions capable of handling the complexities inherent in stochastic differential and partial differential equations. Also, numerical challenges associated with these methods are discussed and the solution strategies are also suggested.

1. Introduction

Mathematical modelling of complex system is a sophisticated way to understand and simulate intricate relationships in various fields including engineering, science, economics, finance and biology [1, 2, 3, 4]. Stochastic Differential Equations (SDEs) and Stochastic Partial Differential Equations (SPDEs) play critical role in understanding the behaviour of these systems. These systems having dynamic interactions and a large number of variables, necessitating advanced mathematical and computational models for proper representation. To capture the complexities and uncertainties inherent in these systems, researchers employ mathematical frameworks using SDEs and SPDEs [5, 6].

Many physical processes are simulated using mathematical models that lead to partial differential deterministic equations. Due to lack of data on parameters and initial data, a system’s behaviour may differ from the ideal deterministic description. Adding random inputs, such as random variables or stochastic processes, to make a system’s explanation more realistic is one way to deal with this issue and compensate for a lack of information. As a result of this, SDEs and SPDEs are formulated.

In financial mathematics, these equations are a central part of modelling; these are used to model the key features such as interest rate, asset prices, risks and their volatility [7, 8, 9]. In applied mathematics, these equations are useful for modelling transport equations [10], quantum mechanics [11], wave propagation in random media [12], neuroscience [13], population genetics and mathematical biology [14]. These equations are also being used in pure mathematics as a natural object to study stochastic processes and dynamical systems and to develop theoretical frameworks, convergence and stability for the algorithms in the form of lemma and theorems.

Although numerous studies have proved the applicability of solute transport models based on deterministic advection-dispersion equations [15, 16], the success of such deterministic models in describing solute movement in variable field conditions is questionable. The deterministic modelling of solute transport in groundwater contamination involves the simulation of hydrodynamic dispersion, advection and absorption with velocity field. Moreover, these components suffer from significant uncertainty because of scale dependency, even in homogeneous media; their effects are poorer in the case of heterogeneous mediums where micro-variations take place in flow parameters.

These difficulties have been captured and handled in two ways. One way is to increase the dispersion parameter with distance, and the other way is to increase it with time. Another approach to scale dependency and heterogeneity is stochastic modelling in which relevant flow parameters are characterised by a distribution function [17]. Ensemble means of the output variables provide a useful description of the larger-scale spatial and temporal trends in practical applications. Although it is not of great importance to know the exact deviation of local fluctuations from the ensemble mean, it is useful to know their likely range.

Various numerical approaches are used to address the challenges posed by the dynamic and uncertain nature of complex systems. Finite difference approaches, for example, discretize the spatial domain and use difference approximations to obtain solutions [18]. This method works well in situations when the spatial features of a complicated system are carefully considered. Monte Carlo simulations, on the other hand, use random sampling to generate numerical results, making them ideal for capturing the randomness found in many complex systems [19, 20]. Efficiency of numerical methods is often measured in terms of computational resources and accuracy. Some numerical methods, like the Implicit-Explicit (IMEX) schemes, strike a balance between stability and efficiency, making them suitable for solving stiff systems with varying time scales. Karhunen-Loève expansion utilises eigenfunctions to represent stochastic processes. This approach is advantageous in problems with random spatial variations, allowing for a reduction of dimensionality.

Another numerical approach to solve SDEs and SPDEs is to use stochastic spectral methods (SSMs). These methods provide a powerful and flexible way for dealing with complex systems defined by these equations and use spectral representations to efficiently handle the systems’ inherent stochasticity, offering correct solutions while maintaining computational efficiency.

Stochastic spectral methods are intrusive and non-intrusive. Intrusive methods [21, 22] represent the solution as a sum of orthogonal basis functions, and then using these basis functions to approximate the solution to the SDE or SPDE. The coefficients of the expansion are determined by solving a system of equations derived from the original problem, typically using techniques such as Galerkin projection or least-squares regression. One popular intrusive stochastic spectral method is the stochastic finite element method (SFEM). SFEM involves representing the solution as a sum of finite element basis functions, and then using these basis functions to approximate the solution [23, 24].

Non-intrusive spectral methods approximate the solution to the SPDE using a series of basis functions. The coefficients of the basis functions are determined by solving a system of equations derived from the original SPDE. These methods are ‘non-intrusive’ in the sense that they do not require the original deterministic solver to be modified in order to handle stochastic inputs. Instead, they use a surrogate model to represent the solution in terms of a spectral expansion. A popular non-intrusive stochastic spectral method is the polynomial chaos expansion (PCE) method. PCE involves representing the solution as a sum of orthogonal polynomials, such as Legendre or Hermite polynomials. The coefficients of the expansion are determined by projecting the original problem onto the polynomial basis, which can be done using techniques such as Galerkin projection or least-squares regression [25, 26, 27]. An advantage of non-intrusive stochastic spectral methods is their ability to handle high-dimensional problems with a relatively small number of samples. This is because the spectral expansion can capture the global behaviour of the solution, allowing for accurate approximation of the solution using only a few terms.

There is a great need and interest in developing efficient numerical schemes to solve these equations. It is essential and profitable too, especially when dealing with randomness. This chapter explores the existing numerical approaches surrounding SDEs and SPDEs relevant to selected applications. Computational aspects of different stochastic spectral methods are discussed and applied to solve some practical problems such as contaminant transport in porous media. A class of problems of solving SDEs/SPDEs using polynomial chaos expansion methods is discussed. The computational efficiencies of polynomial chaos expansion methods are compared with alternative methods. This chapter also explores Wick calculus, a potent but less commonly used method in the stochastic calculus.

2. Numerical approaches for SDEs and SPDEs

2.1 Direct integration techniques

In stochastic analysis, quantities of interest can be represented by the expected value of certain functionals of the dependent variable, say, y, and denoted as E. The value of this expectation is defined in the form of integrals and the need to compute certain integrals, with respect to the probability measure $P_{\xi }$ as:

$$E\left(f\right)={\int }_{R}f\left(y\left(x\right);x\right)d{P}_{\xi }\left(x\right)={\int }_{R}f\left(y\left(x\right);x\right)p_{\xi }\left(x\right)\mathit{dx}$$E1

where $p\left(\xi \right)$ is the probability density function of $\xi$, and the integration strategies lead to the following estimation:

$$E\left(f\right)\approx Q_i\left(f\right)=\sum _{i=1}^{n}f\left(y\left(x_i\right);x_i\right)w_i$$E2

where the $x_i$ are the points and $w_i\in R$ are the weights of integration. The model response is then evaluated for possible outcomes $\xi =x_i$ of the selected random variables. Finally, it requires the solution of a system having $n$ uncoupled deterministic equations.

$$L\left(y\left(x_i\right);x_i\right)=b\left(x_i\right),i=1,2,\dots \dots \dots n$$E3

2.2 Finite difference methods

Finite difference methods involve discretizing the spatial and temporal domains of the equations. For instance, the Euler-Maruyama method is widely used for SDEs. In this approach, the stochastic process is simulated by iteratively updating the solution at each time step.

Consider the example of modelling stock prices $S_t$ as given in Eq. (4), where the geometric Brownian motion equation can be solved using the Euler-Maruyama method. This numerical scheme is the generalisation of Euler’s method for ODEs to solve SODEs. Consider an Itô SDE [28] of the form,

$$d{S}_t=\mu \left(t,S_t\right)\mathit{dt}+\sigma \left(t,S_t\right)d{W}_t$$E4

where $W_t$ is a standard Wiener process on the interval $\left[0,T\right]$. To solve this equation numerically, we discretise the time interval with an $h$ (a sufficiently small) step size. Symbolically, Eq. (4) can be written as,

$$S_t=S_{t_0}+{\int }_{t_0}^t\mu \left(t,S_t\right)\mathit{dt}+{\int }_{t_0}^t\sigma \left(t,S_t\right)d{W}_t$$E5

and the Euler-Maruyama (E-M) approximation is calculated using the Itô-Taylor expansion.

The first integral on the right-hand side of (Eq. (5)) using an Itô-Taylor expansion is,

$${\int }_{t_0}^t\mu \left(t,S_t\right)\mathit{dt}={\int }_{t_0}^t\mu \left(t_0+\left(t-t_0\right),S_0+S\left(t\right)-S_0\right)\mathit{dt}$$E6

$$=\mathit{h\mu }\left(t_0,S_0\right)+O\left(h^{3/2}\right)$$E7

or simply,

$$S_{n+1}=S_n+\mu \left(t_n,S_n\right){\int }_{t_n}^{t_{n+1}}\mathit{d\tau }+\sigma \left(t_n,S_n\right){\int }_{t_n}^{t_{n+1}}d{W}_{\tau }$$E8

where $n=0,1,\dots \dots \dots ,N-1$ for $t_0<t_1<\dots <t=T$ and $N\in N$. Thus, for an $m$-dimensional SDE, $S^h=\left(S_1^h,S_2^h,\dots ,S_m^h\right)$ and

$$S_{t_{n+1}}^h=S_{t_n}^h+\mathit{h\mu }\left(t_n,S_{t_n}^h\right)+\sigma \left(t_n,S_{t_n}^h\right)∆B_n$$E9

with $∆B_n=B_{t_{n+1}}-B_{t_n},t_n=\mathit{nh}$

$∆B_n$ is the sequence of identically independent Gaussian random variables having a zero mean and $h$ standard deviation.

Eq. (9) can be rewritten as,

$$S_{t_{n+1}}-S_{t_n}=\mathit{h\mu }\left(t_n,S_{t_n}^h\right)+\sqrt{h}\sigma \left(t_n,S_{t_n}^h\right)B_n+O\left(h\right)$$E10

The order condition for h judges the quality of the numerical scheme. Consequently, the variation in step size $h$ defines the convergence rate of the numerical scheme. In numerical approximation, based on model requirements, the two most -commonly used definitions of convergence are strong and weak convergence.

2.2.1 Numerical example

We shall apply the E-M method to a simple population growth model which has an exact solution. The population growth model equation of a species $\left(S\right)$ in time $t$ years is written as:

$$d{S}_t=\alpha S_t\mathit{dt}+\beta S_t{\mathit{dW}}_t$$E11

where $\alpha$ and $\beta$ are the model parameters. The exact solution of this equation, using Itô’s formula, is given by:

$$S_t=S_0{e}^{\left(\left(\alpha -\frac{{\beta }^2}{2}\right)t+\beta W_t\right)}$$E12

Simulation (three runs) of the population growth model for an exact solution, with parameter values $\alpha =0.35$ and $\beta =0.25$, is shown in Figure 1. Initial population is assumed to be $S_0=20$ and $\mathit{\Delta t}=0.01$. The E-M approximation is also plotted for the three runs for comparison in Figure 2 with the same parameter values for all computations.

Figure 1.

Exact solution for population growth model.

Figure 2.

Euler-Maruyama approximation for population growth model.

The simulation results are different for the same parameter values because of the stochastic nature of the model. Thus, it is a good idea to run the simulation numerous times and choose an average value to see how the given system behaves (Figure 3).

Figure 3.

Simulation results for the exact and Euler-Maruyama approximation after a number of runs.

The first graph on the top left corner is the average population growth after five simulations. The results from the E-M scheme differ from the exact solution. However, as the simulation runs are increased, the E-M approximations coincide with the actual solution.

2.3 Monte Carlo (MC) methods

Monte Carlo methods are very popular for handling both SDEs and SPDEs. These methods involve random sampling to approximate solutions. MC methods approximate the outcomes based on random sampling and the law of inertia of large numbers. An example is the Black-Scholes model in finance, where Monte Carlo simulations can estimate option prices by simulating the stochastic evolution of asset prices. Practically, these random samples are pseudo-random samples of the basic random variable $\xi$ and integration error is defined by a Gaussian Random Variable (GRV) (Figure 4).

Figure 4.

Generation of a Gaussian random variate (GRV) for different sample sizes and a comparison with MC integration error and the mean.

A basic MC estimator of an integral of a quantity of interest (expectation), say $\mu$, is given by,

$$\text{Mean}\left(\mu \right)={\int }_{R}f\left(x\right)p\left(x\right)\mathit{dx}=E\left[f\left(x\right)\right]\approx \frac{1}{n}\sum _{i=1}^{n}f\left(x_i\right),$$E13

the variance of the estimator is,

$$var\left(\mu \right)=\frac{var\left[f\left(x\right)\right]}{n},$$

The MC estimator does not work accurately for large variance. Statistical errors in the confidence interval can be too large to believe that the interval is valid. For variance reduction, standard deviation of the estimator equals to $\frac{{\sigma }_t}{\sqrt{n}}$. For small confidence interval, we can reduce the variance. For example, $var\left[f\left(x\right)\right]=1$, and a smaller confidence interval is required, say ${10}^{-2}$;

Thus,

$$var\left[f\left(x\right)\right]={10}^{-4}\left(n\right)$$E14

Thus, we need $n={10}^4$ sample points. Numerous techniques based on the control variate, importance sampling, stratified sampling and anthetic variates are available in literature. Zhang and Karniadakis [29] used different reduction methods to avoid generating a large number of sample points. The convergence rate of the basic MC estimator for a one-dimensional integral is $O\left(\frac{1}{\sqrt{n}}\right)$. While this convergence is free from stochastic dimensions, it is slow. The positive side of this independence confirms the applicability of MC methods in problems with high-stochastic dimensions.

2.4 Quasi monte carlo (QMC) methods

Computing mathematical models using MC methods is expensive due to the requirement of a large number of sample points. One solution is to use variance reduction techniques. An error in these approaches is proportional to $\frac{1}{n}$ as compared to $\frac{1}{\sqrt{n}}$ as in the MC methods. A detailed study of these techniques has been given in [30]. A few low-discrepancy sequences are shown in Figure 5.

Figure 5.

Quasi-random sequences. Quasi-random sequences are used in variance reduction techniques. The word ‘quasi’ signifies that the values of a low-discrepancy sequence is neither random nor pseudo-random, but share some of the characteristics of both of the random variables. This low-discrepancy property is an added benefit of QMC methods. Using a control variate helps to reduce MC variance.

2.5 Stochastic spectral methods

Stochastic spectral methods use spectral representations to efficiently capture the inherent randomness and complexity in systems governed by SDEs and SPDEs. Some of the spectral methods are:

2.5.1 Karhunen-Loève expansion

The Karhunen-Loève expansion is particularly useful for problems with random spatial variations. It uses eigenfunctions to represent stochastic processes and helps in dimensionality reduction. In environmental modelling, this method is being applied to simulate groundwater flow in heterogeneous aquifers, where the spatial variability of hydraulic conductivity is considered as a stochastic field.

2.5.2 Polynomial chaos expansion (PCE) methods

The basic idea of PCE methods is to approximate the solution as a sum of polynomial functions that are orthogonal, with respect to the probability distribution of the random inputs [25]. The coefficients of the expansion are determined by solving a system of equations derived from the original problem, typically using techniques such as Galerkin projection or least-squares regression. A schematic for working process of PCE methods is shown in Figure 6.

Figure 6.

Diagrammatic representation for working process of polynomial chaos expansion methods for SDEs and SPDEs.

The resulting PCE can then be used to estimate the solution at any point in the domain, including those with stochastic inputs. One advantage of PCE methods is in their efficiency to propagate uncertainty through the system, allowing for the efficient estimation of statistics such as mean and variance [31]. PCE methods can also be used to efficiently quantify the effect of input uncertainties on the system response, such as in sensitivity analysis. PCE methods can be used in both intrusive and non-intrusive implementations. While intrusive PCE methods modify the original deterministic solver to handle the stochastic inputs directly, non-intrusive PCE methods use a surrogate model to represent the solution in terms of the PCE.

2.5.3 Generalised polynomial chaos expansion (gPCE) methods

Wiener PCE methods use Hermite polynomials and express the solution of SDE/SPDE in terms of these smoothed orthogonal polynomials. Although these homogeneous Wiener polynomial chaos is quite useful and satisfies certain universal polynomial characteristics, but these polynomials are more useful for Gaussian distributed random input parameters. For non-Gaussian random variables, Hermite polynomial-based chaos expansion is slow. Askey scheme for orthogonal polynomials is given in [32]. It involves mapping the stochastic solution onto a set of polynomials that are orthogonal to the probability distribution of the input parameter.

Thus, instead of finding the random solution of the SDE/SPDE, the problem is reduced to evaluating the expansion coefficients. The basis functions are selected based on the distribution of random parameters. The correspondence between the basis functions and their known distributions is available in the literature [32] and is summarised in Figure 7.

Figure 7.

Equivalence between distributions, basis polynomials and related input random variables.

3. Applications of PCE method

3.1 Development of a meta model for geometric brownian motion (gBM)

A sequence of random variables $W={\left\{W\left(t\right)\right\}}_{t\ge 0}\in R$ in a time domain on a probability space $\left(\Omega ,F,\wp \right)$ is known as a standard Brownian motion.

For a fixed $T\in \left(0,\infty \right)$ and a filteration, let $F_t^W$ be the sigma-algebra generated by these random variables for $0\le t\le T.$ Using Cameron-Martin basis [33] for a countable multi-index set,

$$J=\left\{\left({\alpha }_i^j,i,j\ge 1\right),{\alpha }_i^j\in \left\{0,1,\dots \dots \right\}\right\}\text{with}\alpha !={\prod }_{i,j}{\alpha }_i^j!$$E15

here $i$ and $j$ are the number of Gaussian random variables and the Wiener processes, respectively.

According to the Cameron-Martin theorem, there exists an orthonormal basis $m_i,i\ge 1$ of square integrable space $L^2\left(\left[0,T\right]\right),$ such that a set of independent random variables ${\xi }_{\mathit{ij}}$ is defined as,

$${\xi }_{\mathit{ij}}={\int }_0^T{m}_i\left(\tau \right)d{W}_{\tau }^j$$E16

Now, the set of random variables known as Cameron-Martin basis are defined as,

$$T=\left\{{\xi }_{\alpha },\alpha \in J\right\}$$E17

A space generated by this basis is the Wiener chaos space $w$ and a random variable $f\left(\omega \right)\in w,$ or some realisation $\omega$ is written in terms of Wiener chaos expansion as,

$$f=\sum _{\alpha \in J}{a}_{\alpha }{\xi }_{\alpha }$$E18

With the above setup defined, we now develop a meta model for geometric Brownian motion given by the solution of SDE, defined in Eq. (19).

Again,

$${\mathit{dY}}_t=\mu \left(t,Y_t\right)\mathit{dt}+\sigma \left(t,Y_t\right)d{W}_t$$E19

gBM is the solution of Eq. (19).

Applying the Euler-Maruyama method and using independent increments of the Wiener process and for each increment $m=\left\{0,1,\dots \dots ,M-1\right\}$, we have,

$$Y_{t_{m+1}}=Y_{t_m}$$E20

for $0=t_1<t_2<\dots \dots <t_{M-1}<t_M=T$ and $∆t_m=t_{m+1}-t_{m.}$

For meta modelling, we input this multivariate function in the form,

$$Y_{t_m}=M\left(Y\right)$$E21

$$Y=\left({\xi }_{1,}{\xi }_{2,\dots \dots \dots .},{\xi }_{m-1}\right),{\xi }_{i}N\left(0,1\right)$$

and, for the Euler-Maruyama scheme,

$$W\left(t\right)N\left(0,t\right),t\in \left[0,T\right]$$

This model automatically chooses Hermite polynomials as the basis functions since the random variables are Gaussian. More conveniently, Eq. (20) is written as

$$Y_t=Y_{t_m}=Y_0\prod _{m=0}^{M-1}\left(1+\mu \left(t_m,Y_{t_m}\right)∆t_m+\sigma \left(t_m,Y_{t_m}\right)\sqrt{∆t_m}{\xi }_m\right)$$E22

Using a suitable transformation, gBM converts Eq. (19) into the form of Eq. (21). Using the Itô’s lemma and defining $X_t=\mathit{ln}{S}_t$, we get

$${\mathit{dX}}_t=\left(\mu -\frac{1}{2}\right){\sigma }^2\mathit{dt}+\sigma \left(t,X_t\right)d{W}_t$$E23

Now, for any computed value of $X_t$, we have,

$$S_t=e^{X_t}$$E24

and this output is of the form

$$S_T=M\left(Y,⨀\right)$$E25

where $⨀$ is the transformation.

Using a UQLab framework [34] and feeding input $Y$ as gBM, we develop a meta model. Figure 8 represents the PCE coefficients spectrum based on different methods. The hyperbolic truncation scheme is used with q-norm in the range $0\le q-norm\le 1$. The Sobol sequence sampling strategy is used for all methods. The graphs show that the number of PCE coefficients can be reduced using a suitable evaluation method and the sparse sampling technique. The output of the model is summarised in Table 1.

Figure 8.

Spectrum of the PCE coefficients for the meta modelling of geometric Brownian motion using (a) Gaussian quadrature, (b) least angle regression, (c) subspace, (d) ordinary least-squares regression, (e) orthogonal matching pursuit and (f) Bayesian compressive sensing. NNZ represents the number of nonzero coefficients in the chaos expansion.

Parameters	Method
	Gaussian Quadrature	BCS	LARS	SP	OMP	OLS
Number of input variables	3	3	3	3	3	3
Maximal degree	14	8	6	9	5	5
q-norm	1	1	1	0.75	0.75	0.75
Size of full basis	680	165	84	111	32	32
Size of sparse basis	680	165	84	55	24	20
Full model evaluations	3375	500	150	150	150	150
Quadrature error	2.24E-29	2.54E-19	1.04E-13	8.13E-12	3.52E-07	1.43E-07
Mean value	1.0305	1.0305	1.0305	1.0305	1.0305	1.0305
Standard deviation	0.1554	0.1554	0.1554	0.1554	0.1554	0.1554
Coefficient of variation	15.09%	15.09%	15.09%	15.09%	15.09%	15.09%

Table 1.

Meta model’s input and output parameters.

This model is also validated with a set of 1E4 samples. The results for the validation error using different approaches are provided in Table 2.

Validation parameter	Gaussian Quadrature	BCS	LARS	SP	OMP	OLS
Relative error	4.59E-28	1.56E-18	1.27E-13	1.10E-06	7.54E-11	6.92E-07

Table 2.

Relative error for model validation.

3.2 Stochastic contaminant transport

The stochastic advection dispersion equation for contaminant transport [35] is given below:

$$\frac{\partial C}{\partial t}+\mu .\nabla C=\nabla .\left[\kappa \left(x;\omega \right)\nabla C\right]+F\left(t,x;\omega \right),\left(t,x;\omega \right)\in R^{+}\times D\times \Omega$$E26

with the initial condition $C\left(0,x;\omega \right)=C_0\left(x;\omega \right)$, where $C$ is the concentration of the contaminantwithin the bounded domain $D\in R^d\left(d=1,2\right)$, $R^{+}$ is the time domain $\left(0,T\right]$ and $\Omega$ is the sample space. $F\left(t,x;\omega \right)$ is the random forcing term. We consider $F\left(t,x;\omega \right)=\sigma \frac{\partial W}{\partial t}$, where $\sigma$ is the standard deviation of the stochastic term, and $\frac{\partial W}{\partial t}$ represents the time derivative of a Wiener process which introduces stochastic variability. $\mu$ is the mean advection velocity representing the bulk flow of the fluid; $\kappa \left(x;\omega \right)$ is the dispersion coefficient characterising the rate of diffusion and $C_0\left(x;\omega \right)$ is the initial condition. $\kappa \left(x;\omega \right)$ and $C_0\left(x;\omega \right)$ are both real valued functions in $D$ and $\Omega$ for all random inputs.

In gPCE, random numbers are chosen based on the probability distribution of input variables and corresponding orthogonal basis functions. Choosing $\left\{{\phi }_n\right\},n=0,1,\dots \dots$ as orthogonal polynomials for random variable $\xi$ satisfies the orthogonality condition given. Applying gPCE to Eq. (26), we have

$$\kappa \left(x;\omega \right)=\sum _{i=0}^N{\kappa }_i\left(x,t\right){\phi }_i\left(\xi \right)$$E27

$$C\left(t,x;\xi \right)=\sum _{j=0}^N{C}_j\left(t,x\right){\phi }_j\left(\xi \right)$$E28

and

$$F\left(t,x;\xi \right)=\sum _{j=0}^N{F}_j\left(t,x\right){\phi }_j\left(\xi \right)$$E29

The value of $N$ (the total number of expansion terms) is based on the order of polynomials $\left(n\right)$ used and the dimension of the random space $\left(p\right).$ It is given by,

$$N+1=\frac{\left(n+p\right)!}{n!p!}$$E30

Substituting Eq. (27) to Eq. (29) in the governing Eq. (26), we get

$$\begin{array}{lll}\sum _{j=0}^N& \frac{\partial C_j\left(t,x\right)}{\partial t}& {\phi }_j\left(\xi \right)+\mu \sum _{j=0}^N\nabla C_j\left(t,x\right){\phi }_j\left(\xi \right)\\ & & =\nabla .\left[\sum _{i=0}^N{\kappa }_i\left(x\right){\phi }_i\left(\xi \right)\nabla \left(\sum _{j=0}^N{C}_j\left(t,x\right){\phi }_j\left(\xi \right)\right)\right]+\sum _{j=0}^N{F}_j\left(t,x\right){\phi }_j\left(\xi \right)\end{array}$$E31

Rewriting the above equation,

$$\begin{array}{lll}\sum _{j=0}^N& \frac{\partial C_j\left(t,x\right)}{\partial t}& {\phi }_j+\mu \sum _{j=0}^N\nabla C_j\left(t,x\right){\phi }_j=\left[\sum _{i=0}^N\sum _{j=0}^N\nabla .\left[{\kappa }_i\left(x\right)\nabla C_j\left(t,x\right)\right]\right]{\phi }_i{\phi }_j\\ & & +\sum _{j=0}^N{F}_j\left(t,x\right){\phi }_j\end{array}$$E32

Next, applying the Galerkin projection onto each polynomial basis ensures that the induced error is orthogonal to functional space spanned by truncated finite dimensional basis functions ${\phi }_n$. Projecting ${\phi }_m$ for each $m=\left\{0,1,\dots \dots ,N\right\}$ and applying the orthogonality conditions, we get

$$\frac{\partial C_m}{\partial t}+\mu \nabla C_m=\sum _{i=0}^N\sum _{j=0}^N\nabla .\left[{\kappa }_i\left(x\right)\nabla C_j\right]E_{\mathit{imj}}+F_m$$E33

where E_imj$=<{\phi }_i{\phi }_k{\phi }_j>$ and is calculated as $E_{\mathit{imj}}=\frac{<{\phi }_i{\phi }_j,{\phi }_m>}{<{\phi }_m,{\phi }_m>}$ by setting up the definition of ${\phi }_n$. The above equation can be written as,

$$\frac{\partial C_m}{\partial t}+\mu \nabla C_m=\sum _{j=0}^N\nabla .\left[a_{\mathit{mj}}\left(x\right)\nabla C_j\right]+F_m$$E34

or

$$\frac{\partial C_m}{\partial t}=\sum _{j=0}^N\left[a_{\mathit{mj}}\left(x\right){\nabla }^2{C}_j+b_{\mathit{mj}}\left(x\right).\nabla C_m\right]-\mu .\nabla C_m+F_m,m=0,1,\dots \dots \dots ,N$$E35

where,

$$a_{\mathit{mj}}\left(x\right)=\sum _{i=0}^N{\kappa }_i\left(x\right){{\rm E}}_{\mathit{imj}}$$E36

and,

$$b_{\mathit{mj}}\left(x\right)=\sum _{i=0}^N\nabla {\kappa }_i\left(x\right){{\rm E}}_{\mathit{imj}}=\nabla a_{\mathit{mj}}\left(x\right)$$E37

Equation (35) is a set of (N + 1) coupled deterministic equations with an initial condition,

$$C\left(0,x;\omega \right)=\sum _{j=0}^N{C}_j\left(0,x\right)\phi \left(\xi \right)=C_0\left(x;\omega \right)$$

where the coefficients $C_j\left(0,x\right)=C_0\left(x\right)$ in the expansion are deterministic and we can obtain the initial and boundary conditions for each expanded (Eq. (26)). For numerical simulation, we considered a one-dimensional case and the equation,

$$\frac{\mathit{dC}\left(t,x;\omega \right)}{\mathit{dt}}+\mu \frac{\mathit{dC}\left(t,x;\omega \right)}{\mathit{dx}}=\frac{d}{\mathit{dt}}\left[\kappa \left(x;\omega \right)\frac{\mathit{dC}\left(t,x;\omega \right)}{\mathit{dx}}\right]+F\left(t,x;\omega \right)$$E38

with the deterministic initial condition $C\left(0,x;\omega \right)=1$ and boundary conditions $x\in \left[0,1\right)$. Without a loss of generality, randomness in transport is characterised by considering the random forcing term as $F\left(t,x;\omega \right)=\sigma \frac{\partial W}{\partial t}$ with the Wiener process term and $\kappa \left(x;\omega \right)$ is characterised with $\kappa \left(x;\omega \right)=\xi$. While applying MC simulation, 10⁵ random numbers are generated. For polynomial chaos expansion, the randomness in $\sigma$ is captured through the polynomial basis functions and the Wiener process. Figures 9 and 10 show the concentration profile using the MC method with 10⁴ samples and the PCE method of polynomial order 3 for $C\left(0,x;\omega \right)=1.$

Figure 9.

Simulation of concentration using MC method.

Figure 10.

Simulation of concentration using PCE method.

Time domain $R^{+}$ is (0,1] with time step $∆t=0.01$, $C\left(0,x;\omega \right)=1,$ space domain $x\in \left[0,1\right)$ and number of random samples $10000$. Hermite polynomials are of order 3 and maximum order of polynomial chaos expansion is $5$. The mean and standard deviation of the Gaussian distribution are 0 and 1, respectively. Figure 10 shows that PCE method is more efficient in capturing the randomness in $C\left(t,x;\omega \right)$ compared to MC simulation.

3.2.1 Effect of boundary conditions on solution

The choice of boundary conditions plays an important role in determining the appropriate order of polynomial chaos expansion (PCE) for computer simulations.

Graphs for the Gaussian initial condition $C\left(0,x;\omega \right)=\frac{1}{{\sigma }_G\sqrt{2\pi }}\mathit{exp}\left({\left(x-\mu \right)}^2,/2{\sigma }_G\right)$ for different order of Hermite polynomials are shown in Figure 11(a)–(d). Boundary conditions determine the behaviour of a random process at the boundaries of the domain of interest. For example, imposing periodic or reflection boundary conditions can have a significant impact on the statistical properties of the random field. In cases where the boundary conditions introduce additional sources of randomness or nonlinearity, higher-order PCE may be needed to accurately capture these effects.

Figure 11.

Simulation of concentration $C$ using the PCE method with the Gaussian initial condition and various order of polynomial chaos expansion. Time domain $R^{+}$ is (0,1], time step $∆t=\mathit{0.01}$, Hermite polynomials of order 3. The mean and standard deviation of Gaussian distribution is $\mu =\mathit{0.5}$ and standard deviation ${\sigma }_G=\mathit{0.1}$.

In contrast, simpler boundary conditions, such as Dirichlet or Neumann conditions, can lead to smoother responses, allowing lower-order PCEs to provide accurate approximations. Therefore, choosing the appropriate PCE order is closely related to understanding the boundary conditions and their influence on the stochastic behaviour of the system, thereby ensuring that the computational model faithfully represents the physical reality of the problem. The graphs in Figure 11 show that lower-order polynomials are capable of capturing the Gaussian distribution more accurately than higher-order polynomials.

4. Wick product

The Wick product was initially used as a renormalisation process. T. Hida and N. Ikeda [36] proposed the Wick product for stochastic analysis. Later on, this concept was expanded [37] to include white noise in the Wick products of stochastic distributions and applied Gaussian Wick calculus. A detailed and clear overview of the Wick product and Wick renormalisation is explained in [38].

4.1 Motivation for Wick multiplication (renormalisation point of view)

White noise $W\left(t\right)$ can be represented as an element of Hida dual space ${\left(L^2\right)}^{-1}$ as,

$$W\left(t\right)=\underset{\mathit{\Delta t}\to 0}{\mathit{lim}}\frac{1}{\mathit{\Delta t}}\underset{t}{\overset{t+\mathit{\Delta t}}{\int }}d{B}_u$$E39

$${\int }_R{\delta }_{t}d{B}_u$$E40

where ${\delta }_t$ is the Dirac measure of $t$ that belongs to the dual Sobolov space. Using Itô’s formula,

$${\left(W_t\right)}^2\approx {\left(\frac{1}{\mathit{\Delta t}}\underset{t}{\overset{t+\mathit{\Delta t}}{\int }}d{B}_u\right)}^2$$E41

$$\frac{2}{{\left(\mathit{\Delta t}\right)}^2}\underset{t}{\overset{t+\mathit{\Delta t}}{\int }}\left(\left(\underset{t}{\overset{v}{\int }}d{B}_u\right)d{B}_v+\frac{1}{\mathit{\Delta t}}\right)$$E42

The additive re-normalisation of the above term is,

$$\frac{2}{{\left(\mathit{\Delta t}\right)}^2}\underset{t}{\overset{t+\mathit{\Delta t}}{\int }}\left(\left(\underset{t}{\overset{v}{\int }}d{B}_u\right)d{B}_v\right)$$

This motivates the definition of,

$$\int {\delta }_t\mathit{dB}◊{\delta }_t\mathit{dB}=\int {\delta }_t\otimes {\delta }_t{\mathit{dB}}^{\otimes 2}$$E43

Here, $\otimes$ is the symmetrized tensor product. The generalised structure of the solutions is transferred to the stochastic component, and the SPDEs are interpreted in the typical strong sense with respect to the time parameter $t$ and the spatial parameter $x\in R^d$. An important aspect of using the Wick product is the ability to define a unique white noise process as a mathematical object, such as the time derivative of the Brownian motion. SDEs can therefore be solved as true differential equations, which is how they are typically understood, rather than just as integral equations.

The Wick product is used in the context of white noise analysis to solve the issue of pointwise multiplication of generalised functions. It is specified in the test and generalised stochastic function of Kondratiev spaces. It is also useful to replace all products with Wick products and all nonlinear functions with their Wick equivalents in order to solve the multiplicative noise and nonlinear equations because the Wick product is a product in the algebraic sense. Although the Wick product performs very smoothly as a mathematical object, care should be taken when using it.

4.2 The Wick product in stochastic analysis

If ${\left(S\right)}^{*}$ is a Hida distribution space and $f\left(\omega \right)$ and $g\left(\omega \right)$ are two elements of ${\left(S\right)}^{*}$ defined as,

$$f\left(\omega \right)={\sum }_{\alpha }{a}_{\alpha }{H}_{\alpha }$$E44

and

$$g\left(\omega \right)={\sum }_{\beta }{b}_{\beta }{H}_{\beta }$$E45

where $\alpha$ and $\beta$ are the multi-indices, then the Wick product of $f\left(\omega \right)$ and $g\left(w\right)$ is defined as,

$$f\left(\omega \right)◊g\left(\omega \right)={\sum }_{\alpha ,\beta }{a}_{\alpha }{b}_{\beta }{H}_{\alpha +\beta }$$E46

Also, $f\left(\omega \right)◊g\left(\omega \right)\in {\left(S\right)}^{*}$. A detailed explanation and application of the Wick product can be found in the extant literature [39].

4.3 Properties of the Wick product

The following are some of the properties of the Wick product used in stochastic analysis:

$$W^{◊2}\left(t\right)=W^2\left(t\right)-t,t\ge 0$$

1. For a random process $X$, $E\left[{exp}^{◊},\left(X\right)\right]=\text{expE}\left[\mathrm{X}\right]$

   $$\mathit{E}\left[\mathit{X}◊\mathit{Y}\right]=\mathit{E}\left[\mathit{X}\right]\mathit{E}\left[\mathit{Y}\right]$$

Chain rule: for an analytic function $g:C\to C,i.e.,g\left(R\right)\subseteq R,$

$$\frac{d}{\mathit{dt}}\left[g^{◊}\left(x\left(t\right)\right)\right]={(g’)}^{◊}x{\left(\right(t)}^{◊}x’\left(t\right))$$E47

also,

$$\left({\int }_{R}g\left(x\right)\mathit{dW}\left(x\right)\right)◊\left({\int }_{R}h\left(y\right)\mathit{dW}\left(y\right)\right)=\left({\int }_{R}g\left(x\right)\mathit{dW}\left(x\right)\right)\left({\int }_{R}h\left(y\right)\mathit{dW}\left(y\right)\right)-{\int }_{R}g\left(t\right)h\left(t\right)\mathit{dt},$$E48

$\forall g,h\in L^2\left(R\right)$ and are deterministic.

1. In general, the $n$th Wick exponential of $\vartheta$ is defined as

$${\vartheta }^{◊n}={||f||}^n{H}_n\left(\frac{\vartheta }{||f||}\right)$$E49

where $\vartheta ={\int }_R\mathit{fdW}$ and hence,

$$W^{◊n}\left(t\right)=t^{\frac{n}{2}}{H}_n\left(\frac{W\left(t\right)}{t^{\frac{1}{2}}}\right),n=0,1,2,\dots \dots \dots$$E50

$${exp}^{◊}\left\{{\int }_R\mathit{fdW}\right\}=exp\left\{{\int }_R\mathit{fdW}-\frac{1}{2}{‖f‖}^2\right\},f\in L^2\left(R\right)$$E51

and,

$${exp}^{◊}W\left(t\right)=exp\left\{W\left(t\right)-\frac{1}{2}t\right\},\forall t\ge 0$$E52

Equations (49) to (52) are very important from computational point of view. These equations help us to understand the basic building block of smoothed white noise process and need of Wick product.

4.4 Computation of wick product

Let us consider (Eq. 52) and compute the Wick product and Wick exponential. First, we will compute least-square norm of a deterministic function $f\left(t\right)$, that is, $f\in L^2\left(R\right)$ given by a second-degree polynomial

$$f\left(t\right)=a+\mathit{bt}+c{t}^2$$E53

$L^2$-norm of $f\left(t\right)$ is

$$‖f‖={\left(\underset{0}{\overset{t}{\int }}{\left(f\left(t\right)\right)}^2\mathit{dt}\right)}^{1/2}$$E54

When $f\left(t\right)=1,‖f‖={\left(\underset{0}{\overset{t}{\int }}\mathit{dt}\right)}^{1/2}=\sqrt{t}$, when $a=0,b=1,c=0;$ $‖f‖={\left(\underset{0}{\overset{t}{\int }}{\left(f\left(t\right)\right)}^2\mathit{dt}\right)}^{1/2}=\sqrt{\frac{t^3}{3}}$ and when $a=1,b=2,c=0.001;$ $‖f‖=\sqrt{2\times {10}^{-7}{t}^5+0.001t^4+1.334t^3+2t^2+t}$ and is plotted in Figure 12.

Figure 12.

Least square norm of a second-degree polynomial function.

Now, let $\vartheta \left(t\right)={\int }_{R}f\left(t\right)\mathit{dW}$; this is an Itô process, and a realisation of $\vartheta$ for the interval $\left[0,1\right]$ is shown in Figure 13.

Figure 13.

Realisation of an Itô process.

Let $\Psi \left(t\right)=\frac{\vartheta }{‖f‖}$ so that using (Eq. 50), we have,

$$\Psi \left(t\right)=\frac{\underset{0}{\overset{t}{\int }}f\left(t\right)\mathit{dW}\left(t\right)}{‖f‖}$$E55

As we have already computed the numerator $\underset{0}{\overset{t}{\int }}f\left(t\right)\mathit{dW}\left(t\right)$ and denominator $‖f‖$, we can compute $\Psi \left(t\right)$, and then, incorporating this $\Psi \left(t\right)$ into (Eq. 50) we can compute ${\vartheta }^{◊n}$ as

$${\vartheta }^{◊n}={‖f‖}^n{H}_n\left(\Psi \left(t\right)\right)$$

Further, using (Eq. 52), we have

$${exp}^{◊}\left\{{\int }_{R}f\left(t\right)\mathit{dW}\left(t\right)\right\}={exp}^{◊}\left\{\vartheta \left(t\right)\right\}=exp\left\{\vartheta \left(t\right)-\frac{1}{2}{‖f‖}^2\right\}$$E56

also, if $f\left(t\right)=1,$

$${exp}^{◊}\left\{{\int }_{R}f\left(t\right)\mathit{dW}\left(t\right)\right\}={exp}^{◊}\left\{W\left(t\right)\right\}=exp\left\{W\left(t\right)-\frac{1}{2}t\right\}$$E57

for $f\left(t\right)=1+2t+0.001t^2$, the graph for Wick exponential (Eq. 55) is plotted in Figure 14.

Figure 14.

Wick exponential of an Itô process $\mathit{\vartheta }\left(\mathit{t}\right)\mathit{=}{\int }_{\mathit{R}}\mathit{f}\left(\mathit{t}\right)\mathit{dW}\left(\mathit{t}\right)$, $\mathit{f}\left(\mathit{t}\right)\mathit{=}\mathit{a}\mathit{+}\mathit{bt}\mathit{+}\mathit{c}{\mathit{t}}^{\mathit{2}}$ with $\mathit{a}\mathit{=}\mathit{1}\mathit{,}\mathit{b}\mathit{=}\mathit{2}\mathit{,}$ and $\mathit{c}\mathit{=}\mathit{0.001}$.

5. Discussion and conclusion

Numerical approaches offer a systematic and effective way to simulate the dynamic behaviour of complex systems, allowing academics and practitioners to obtain insights into their complex dynamics. Numerical methods are vital tools for modelling and analysis of these systems, where the interaction of many elements creates stochasticity. Finite difference/element methods, Monte Carlo simulations and spectral methods are all adaptable approaches to dealing with SDEs and SPDEs. These methods assist the investigation of complex events in a variety of domains by providing computationally feasible solutions.

The ability of numerical approaches to deal with high-dimensional problems and spatial changes is very important when studying complex systems. For example, Stochastic Collocation and Galerkin approaches enable researchers to address uncertainty in several dimensions and spatial domains. This adaptability is essential for reflecting the complexities of real-world systems where uncertainties emerge in a variety of ways. Furthermore, as computer capabilities improve, numerical approaches play an increasingly important role in enhancing our knowledge of complex systems. These give a practical way to investigate and simulate stochastic systems, allowing for the exploration of problems that would otherwise be intractable using analytical methods alone.

The toolbox of stochastic spectral approaches is extensive and diverse. While stochastic spectral approaches provide an impressive set of tools for dealing with uncertainty, these approaches are not without limitations. One of the most significant issues is the computational intensity, particularly when dealing with high-dimensional situations. Monte Carlo simulations, for example, may require a large number of samples to reach convergence, making them computationally expensive. Similarly, polynomial chaos expansions have difficulties in adequately capturing highly nonlinear dynamics, frequently demanding an increased number of terms or advanced techniques to accurately approximate the solution. As a result, these methods may struggle to give real-time solutions for complicated systems, limiting their application in situations requiring quick decision-making.

Another drawback comes from the inbuilt assumptions of the modelling process. Numerical methods often assume certain statistical properties of the underlying uncertainties, such as square integrable spaces, specific probability distributions or correlations. These assumptions, while simplifying the problem, may not always align with the true nature of the system, leading to model inaccuracies. Also, the curse of dimensionality remains a serious issue, especially in high-dimensional systems, where it can exponentially increase the computational cost. As a result, when applied to real-world scenarios with large parameter spaces, numerical methods may find practical challenges. To address these constraints, continued research into more efficient algorithms, adaptive strategies and novel methodologies to improve the adaptability and dependability of numerical approaches in dealing with complex, high-dimensional and nonlinear systems is required.