Solution of the Stochastic Heat Equation with Nonlinear Losses Using Wiener-Hermite Expansion

Abstract

In the current work, the Wiener-Hermite expansion (WHE) is used to solve the stochastic heat equation with nonlinear losses. WHE is used to deduce the equivalent deterministic system up to third order accuracy. The solution of the equivalent deterministic system is obtained using different techniques numerically and analytically. The finite-volume method (FVM) with Pickard iteration is used to solve the equivalent system iteratively. The WHE with perturbation technique (WHEP) is applied to deduce more simple and decoupled equivalent deterministic system that can be solved numerically without iterations. The system resulting from WHEP technique is solved also analytically using the eigenfunction expansion technique. The Monte-Carlo simulations (MCS) are performed to get the statistical properties of the stochastic solution and to verify other solution techniques. The results show that higher-order solutions are essential especially in case of nonlinearities where non-Gaussian effects cannot be neglected. The comparisons show the efficiency of the numerical WHE and WHEP techniques in solving stochastic nonlinear PDEs compared with the analytical solution and MCS.

1. Introduction

The heat equation is of fundamental importance in diverse scientific fields. In mathematics, it is the prototypical parabolic partial differential equation. In probability theory, the heat equation is connected with the study of Brownian motion via the Fokker-Planck equation. In financial mathematics it is used to solve the Black-Scholes partial differential equation. The diffusion equation, a more general version of the heat equation, arises in connection with the study of chemical diffusion and other related processes [1].

The existence and uniqueness of stochastic PDEs (SPDEs) were studied extensively in the literature; see, for example, [2–6]. In [4, 6] the Green function is used to get a solution of the one-dimensional heat equation with white-noise source. The white noise was approximated and the authors showed that the solutions converge weakly to the solution of the original equation with white noise. Also in [6], the authors studied the existence, uniqueness, and continuity of the nonlinear diffusion equation in case the nonlinearity term is of polynomial type. In [7] the long term behavior of the solution of the stochastic heat equation with spatially colored random force is considered. The conditions under which the stochastic PDE admits a unique solution were outlined. More generally, the existence of a martingale solution for parabolic stochastic PDEs driven by Poissonian noise is shown in [8].

The solution of SPDEs using WHE has the advantage of converting the problem into a system of equivalent deterministic equations that can be solved efficiently using the standard deterministic methods. The main statistics, such as the mean, covariance, and higher-order statistical moments, can be calculated by simple formulae involving only the deterministic WHE coefficients [9]. The application of the WHE aims at finding a truncated series solution to the solution process of differential equations. The truncated series is composed of two major parts; the first is the Gaussian part which consists of the first two terms, while the rest of the series constitute the non-Gaussian part. In nonlinear cases, difficulties exist in solving the resultant set of deterministic integrodifferential equations obtained from the applications of a set of comprehensive averages on the stochastic integrodifferential equation obtained after the direct application of WHE. Many authors introduced different techniques to solve these difficulties. Among them, the WHEP technique was introduced in [9] using the perturbation technique to simplify and decouple the WHE deterministic system. The WHEP technique was automated and generalized for the case of polynomial type nonlinearity in [10]. The equivalent deterministic system can, using the automation, be obtained for any order and any number of corrections. This enhancement allows for deducing and hence getting higher-order approximations.

The WHEP technique was used to solve different stochastic linear and nonlinear PDEs; see, for example, [11–14]. The one-dimensional diffusion equation was solved numerically in [13] using both WHE and WHEP techniques. It was found from the comparisons there that the WHE solution is the limit of the WHEP solution as the number of corrections approaches infinity. The numerical WHE solution can be integrated with Grunwald-Letnikov numerical approximation of the fractional derivatives to solve stochastic PDEs with fractional or variable order derivatives [12].

Other techniques can be used in solving the stochastic PDEs with nonlinearities, such as the homotopy perturbation (HPM) technique, homotopy analysis (HAM) technique, and the Pickard’s iteration technique. Comparisons between these techniques with the WHEP technique are done in [14]. The WHEP technique was found to be the most efficient technique among other techniques. The stochastic heat equation with nonlinear losses was solved with the HAM technique in [15] and is compared with the WHEP solutions. It was found also that the WHEP technique is simpler and more efficient than HAM technique.

The stochastic heat equation is solved using WHEP technique [11, 14, 15] up to first order only. The first order solutions are suitable for linear problems but in case of nonlinear problems, higher-order solutions are essential. The non-Gaussian behavior of the solution can be detected only with higher-order terms. Higher-order solutions are obtained using automated WHEP technique for the oscillatory equation in [16].

The new contribution in this work is the higher-order solutions of the stochastic heat equation with nonlinear losses of polynomial type. The non-Gaussian part of the solution is evaluated compared with the total solution to highlight the effect of nonlinearities that cannot be neglected. Additionally, the solution is verified using four different techniques: the WHE with numerical Pickard’s iterations, WHEP with numerical estimation, analytical eigenfunction expansion solved using Mathematica, and the Monte-Carlo simulations (MCS).

This paper is organized as follows: the problem is introduced in Section 2. The WHE technique is explained in Section 3 and the equivalent deterministic system is deduced. The numerical Pickard’s iterations are explained also in this section. In Section 4, the perturbation technique is added to the WHE to deduce an extended set of deterministic equations. Section 5 describes the analytical solutions and how to evaluate it using Mathematica. Section 6 describes the MCS of the original stochastic heat equation. Section 7 shows the results and comparisons of the different techniques.

2. Problem Formulation

Consider the following heat equation with nonlinear losses and zero boundary conditions in the form 𝜕𝑢𝜕𝜕𝑡−𝛼²𝑢𝜕𝑥²=−𝜆𝑢^𝑛(𝑡)𝑁();()∈()×(),𝑢()()𝑢()(𝑥),+𝑓𝑥;𝑤𝑡,𝑥0,∞0,𝑡,0=𝑢𝑡,=0,0,𝑥=𝜑 (1) where 𝛼 is the thermal diffusivity and 𝜔 is the random outcome of a triple probability space (Ω,𝜒,𝑃) in which Ω is a sample space, 𝜒 is 𝜎-field associated with Ω, and 𝑃 is a probability measure. The nonlinear losses term 𝑢^𝑛 is scaled (strengthened) by a deterministic scalar 𝜆. The term 𝑁(𝑥;𝑤) is a space white-noise excitation term that has the property 𝐸[𝑁(𝑥;𝑤),𝑁(𝑥₁;𝑤)]=𝛿(𝑥−𝑥₁). The function 𝑓(𝑡) is a deterministic envelope function. The initial condition 𝜑(𝑥) will be assumed to be deterministic. The nonlinear losses term 𝑢^𝑛 may account for the radiative heat loss which depends on the excess heat 𝑢 compared with the surroundings. At low excess temperatures, the radiative loss may be approximated linearly (proportional to 𝑢) while at high excess temperature, the Stefan-Boltzmann law gives a net radiative heat loss proportional to 𝑢⁴.

According to theorem  4.1 in [17], (1) has a unique finite mean solution if the functions 𝜆𝑢^𝑛 and 𝑓(𝑡) are globally Lipschitz and if the initial condition 𝜑(𝑥)∈𝐶([0,]). In [18], the authors provided the sufficient conditions to ensure that the real valued mild solution of the SPDE (1) after approximating the white-noise term converges in law to the solution of the original white-noise driven equation in (1).

Define the linear operator 𝐿 as 𝜕𝐿=𝜕𝜕𝑡−𝛼²𝜕𝑥². (2) In the current work and without loss of generality, a unit thermal diffusivity (𝛼=1) and second degree nonlinear losses (𝑛=2) are assumed. Now, (1) can be rewritten as 𝐿[𝑢(𝑡,𝑥;𝑤)]=−𝜆𝑢²+𝑓(𝑡)𝑁(𝑥;𝑤), (3) with the initial and boundary conditions described in (1).

3. Description of the WHE Technique

In the WHE technique, the stochastic response function 𝑢(𝑡,𝑥;𝑤) is expanded as [10] 𝑢(𝑡,𝑥;𝑤)=^∞_{𝑘=0𝑅𝑘}𝑢^(𝑘)𝑡,𝑥;𝑥₁,𝑥₂,…,𝑥_𝑘×𝐻^(𝑘)𝑥₁,𝑥₂,…,𝑥_𝑘𝑑𝜏_𝑘, (4) where 𝑢^(𝑘)(𝑡,𝑥;𝑥₁,𝑥₂,…,𝑥_𝑘) is the 𝑘th deterministic kernel of 𝑢(𝑡,𝑥;𝑤), 𝑑𝜏_𝑘=𝑑𝑥₁𝑑𝑥₂⋯𝑑𝑥_𝑘, and ∫_𝑅𝑘 is a 𝑘-dimensional integral over the disposable variables 𝑥₁,𝑥₂,…,𝑥_𝑘. The functional 𝐻^(𝑘)(𝑥₁,𝑥₂,…,𝑥_𝑘) is the 𝑘th order Wiener-Hermite functional. The WHE technique for general nonlinear exponent (𝑛) and general order (𝑚) can be summarized in the following steps [10]. For the sake of brevity, the space and time variables (𝑡,𝑥) and 𝜔 will be omitted and only disposable variables 𝑥₁,𝑥₂,𝑥₃,… are written: (1)truncate the expansion (4) to contain only 𝑚+1 kernels 𝑢^(𝑗);0≤𝑗≤𝑚; that is, 𝑢(𝑡,𝑥;𝑤)=𝑢⁽⁰⁾(𝑡)+∑^𝑚_𝑗=1∫_𝑅𝑗𝑢^(𝑗)𝐻^(𝑗)𝑑𝜏_𝑗,(2)substitute the truncated expansion into  (3),(3)use the multinomial theorem to expand the nonlinear term 𝑢^𝑛;𝑛=2,(4)multiply by 𝐻^(𝑗);0≤𝑗≤𝑚 and then apply the ensemble average.This will result in (𝑚+1) equations in the deterministic kernels 𝑢^(𝑗);0≤𝑗≤𝑚 as 𝐿𝑢^(𝑗)𝜆=−𝑗!_𝑔𝑐_𝑔𝑅𝑧^𝑚_𝑖=0𝑢^{(𝑖)𝑘𝑖𝑔}𝐸^𝑗_𝑔𝑑𝜏_𝑧+𝛿_𝑗1𝑓(𝑡)𝛿𝑗!𝑥−𝑥₁;0≤𝑗≤𝑚. (5) The expectations 𝐸^𝑗_𝑔 are computed as 𝐸^𝑗_𝑔=⟨𝐻^(𝑗)∏^𝑚_𝑖=0(𝐻^(𝑖))^𝑘𝑗𝑔⟩. It was explained in [10] how to get 𝐸^𝑗_𝑔 in terms of the Dirac delta functions and then use them to reduce the order of integration. The Kronecker delta function 𝛿_𝑗1 equals one when 𝑗=1 and zero otherwise. The counter 𝑔 in the summation in the right-hand side of (5) runs over all the (_𝑛^𝑛+𝑚) combinations of the positive integers 𝑘⁰_𝑔,𝑘¹_𝑔,…,𝑘^𝑚_𝑔 such that ∑^𝑚_𝑖=0𝑘^𝑖_𝑔=𝑛. In case of 𝑚=1, the Gaussian part of the solution is obtained. Larger values of 𝑚 will account for the non-Gaussian part of the solution.

After solving (5) for the kernels 𝑢^(𝑗);0≤𝑗≤𝑚, the expectation and variance of the solution are obtained as [10] 𝐸[𝑢(𝑡,𝑥;𝑤)]=𝑢⁽⁰⁾,Var[𝑢(𝑡,𝑥;𝑤)]=^𝑚_𝑘=1(𝑘!)_𝑅𝑘𝑢^(𝑘)2𝑑𝜏_𝑘. (6) In case of quadratic nonlinearity (𝑛=2) and 3rd order approximation (𝑚=3), (5) will expand to the following deterministic system: 𝐿𝑢⁽⁰⁾=−𝜆𝑢⁽⁰⁾²−𝜆_𝑅𝑢⁽¹⁾𝑥₁²𝑑𝑥₁−2𝜆_𝑅2𝑢⁽²⁾𝑥₁,𝑥₂²𝑑𝑥₁𝑑𝑥₂−6𝜆_𝑅3𝑢⁽³⁾𝑥₁,𝑥₂,𝑥₃²𝑑𝑥₁𝑑𝑥₂𝑑𝑥₃,𝐿𝑢⁽¹⁾=𝑓(𝑡)𝛿𝑥−𝑥₁−2𝜆𝑢⁽⁰⁾𝑢⁽¹⁾𝑥₁−4𝜆𝑢⁽¹⁾𝑥₂𝑢⁽²⁾𝑥₁,𝑥₂𝑑𝑥₂−12𝜆_𝑅2𝑢⁽²⁾𝑥₂,𝑥₃𝑢⁽³⁾𝑥₁,𝑥₂,𝑥₃𝑑𝑥₂𝑑𝑥₃,𝐿𝑢⁽²⁾=−2𝜆𝑢⁽⁰⁾𝑢⁽²⁾𝑥₁,𝑥₂−𝜆𝑢⁽¹⁾𝑥₁𝑢⁽¹⁾𝑥₂−6𝜆_𝑅𝑢⁽¹⁾𝑥₃𝑢⁽³⁾𝑥₁,𝑥₂,𝑥₃𝑑𝑥₃−4𝜆_𝑅𝑢⁽²⁾𝑥₁,𝑥₃𝑢⁽²⁾𝑥₂,𝑥₃𝑑𝑥₃−18𝜆_𝑅2𝑢⁽³⁾𝑥₁,𝑥₃,𝑥₄𝑢⁽³⁾𝑥₂,𝑥₃,𝑥₄𝑑𝑥₃𝑑𝑥₄,𝐿𝑢⁽³⁾=−2𝜆𝑢⁽⁰⁾𝑢⁽³⁾𝑥₁,𝑥₂,𝑥₃−2𝜆𝑢⁽¹⁾𝑥₁𝑢⁽²⁾𝑥₂,𝑥₃−4𝜆_𝑅𝑢⁽²⁾𝑥₁,𝑥₄𝑢⁽³⁾𝑥₂,𝑥₃,𝑥₄𝑑𝑥₄−4𝜆_𝑅𝑢⁽²⁾𝑥₂,𝑥₄𝑢⁽³⁾𝑥₁,𝑥₃,𝑥₄𝑑𝑥₄−4𝜆_𝑅𝑢⁽²⁾𝑥₃,𝑥₄𝑢⁽³⁾𝑥₁,𝑥₂,𝑥₄𝑑𝑥₄. (7) The zero boundary condition is applied to all kernels in (7). The initial condition in (1) will be applied only to 𝑢⁽⁰⁾ while all other kernels will have zero initial conditions.

The equivalent deterministic system (7) is coupled and nonlinear. The analytic solution of (7) is not easy to obtain even with the available computer packages. The solution can be obtained numerically using suitable numerical technique. The FDM is not suitable in this case as there is a term with Dirac delta function that requires an integral scheme such as FEM or FVM. The nonlinearity can be manipulated with Pickard’s iterations. The nonlinear terms are written in the right-hand side and computed explicitly (from the previous time iteration). The equivalent deterministic system in (7) can be written in the following model form for any 𝑣=𝑢^(𝑗);0≤𝑗≤𝑚: 𝐿(𝑣)=𝐹𝑣,𝑢⁽⁰⁾,…,𝑢^(𝑚);𝑣0,𝑥;𝑥₁,…,𝑥_𝑗=𝑣₀. (8) Pickard’s iteration takes the form 𝑣^𝑘+1=𝐿⁻¹𝐹𝑣,𝑢⁽⁰⁾,…,𝑢^(𝑚)𝑘, (9) where 𝑘 is the time level. Other fixed-point iterations can be also used, for example, Krasnoselski-Schaefer, Mann, and Ishikawa [19], but in our cases, Pickard’s iteration was found to converge faster than other iterations. This is due to the contractive properties of the operator 𝐿⁻¹ in (9). If the operator 𝐿⁻¹ is not strictly contractive, the other fixed-point techniques should be used because Pickard’s iterations will not converge to a fixed point in this case.

There is a classical general existence theory of fixed points for mappings satisfying compactness conditions associated with the names of Brower, Schander, Leray, and so forth, as well as an abundant literature of metrical fixed point theorems, which establish the existence (and uniqueness) of fixed points for maps satisfying a variety of contractive conditions [19]. The first basic result is the classical Picard-Banach-Caccioppoli principle which guarantees the convergence of (9) to a fixed point if the right-hand side operator 𝐿⁻¹ is a strict contraction or satisfies a certain generalized strong contraction condition. In our case, the contraction condition of 𝐿⁻¹ depends on the used numerical scheme and the value of the nonlinearity strength (𝜆).

We can apply the FVM with the semi-implicit Crank-Nicolson scheme. Crank-Nicolson scheme has the advantage that it is unconditionally stable [20], assuming that Δ𝑡 and Δ𝑥 are the widths of the subdivisions in the time and space directions, respectively. Applying this numerical scheme to (9) will result in the following tridiagonal system: −0.5𝑟𝑣^𝑘+1_𝑖−1+(1+𝑟)𝑣_𝑖^𝑘+1−0.5𝑟𝑣^𝑘+1_𝑖+1=0.5𝑟𝑣^𝑘_𝑖−1+(1−𝑟)𝑣^𝑘_𝑖−0.5𝑟𝑣^𝑘_𝑖+1+Δ𝑡Δ𝑥𝐹^𝑘_𝑖, (10) where 𝑣^𝑘_𝑖 is the kernel 𝑣 at time level 𝑘 and space subdivision 𝑖. The factor 𝑟 is defined as 𝑟=Δ𝑡/(Δ𝑥)², and 𝐹^𝑘_𝑖 is defined as 𝐹^𝑘_𝑖=^𝑥𝑖_𝑥^+Δ𝑥/2_{𝑖−Δ𝑥/2}𝐹^𝑘𝑣,𝑢⁽⁰⁾,…,𝑢^(𝑚)𝑑𝑥. (11) The above described scheme in (10) has a truncation error of 𝑂(Δ𝑡²,Δ𝑥²) [20]. The computations are to be repeated until the residual becomes below a specified stopping criterion. In the current work, the residual is computed as Residual=^𝑚_𝑗=0𝑢^(𝑗)_𝑘+1−𝑢_𝑘^(𝑗). (12)

4. The WHEP Technique

The perturbation technique can be applied to the coupled integrodifferential system in (7) to produce simpler decoupled deterministic system [10]. It can be applied up to a certain number of corrections 𝑁𝐶. In the WHEP technique, the following two steps will be added to the above 4 steps of the WHE technique.(5)For each kernel 𝑢^(𝑗);0≤𝑗≤𝑚, apply the perturbation technique up to 𝑁𝐶 corrections; that is, substitute 𝑢^(𝑗)=∑^𝑁𝐶_𝑖=0𝜆^𝑖𝑢_𝑖^(𝑗) in (7).(6)Compare the coefficients of 𝜆^𝑘;0≤𝑘≤𝑁𝐶 in both sides to get 𝑁𝐶+1 equations for each kernel 𝑢^(𝑗);0≤𝑗≤𝑚.This will result in (𝑚+1)(𝑁𝐶+1) deterministic equations. The WHEP technique has the advantage of decoupling the equivalent deterministic system resulting from WHE technique. This means that the numerical solution can be applied directly without iterations.

Applying the WHEP technique to the system in (7), the final set of equivalent deterministic equations for 𝑁𝐶=2 will be 𝐿𝑢₀⁽⁰⁾𝐿𝑢=0,₀₍₁₎=𝑓(𝑡)𝛿𝑥−𝑥₁,𝐿𝑢₀⁽²⁾𝐿𝑢=0,₀₍₃₎𝐿𝑢=0,₁₍₀₎=−𝑢₀⁽⁰⁾²−_𝑅𝑢₀⁽¹⁾𝑥₁²𝑑𝑥₁,𝐿𝑢₁⁽¹⁾=−2𝑢₀⁽⁰⁾𝑢₀⁽¹⁾𝑥₁,𝐿𝑢₁⁽²⁾=−𝑢₀⁽¹⁾𝑥₁𝑢₀⁽¹⁾𝑥₂,𝐿𝑢₁⁽³⁾𝐿𝑢=0,₂₍₀₎=−2𝑢₀⁽⁰⁾𝑢₁⁽⁰⁾−2_𝑅𝑢₀⁽¹⁾𝑥₁𝑢₁⁽¹⁾𝑥₁𝑑𝑥₁,𝐿𝑢₂⁽¹⁾=−2𝑢₀⁽⁰⁾𝑢₁⁽¹⁾𝑥₁−2𝑢₁⁽⁰⁾𝑢₀⁽¹⁾𝑥₁−4_𝑅𝑢₀⁽¹⁾𝑥₂𝑢₁⁽²⁾𝑥₁,𝑥₂𝑑𝑥₂,𝐿𝑢₂⁽²⁾=−2𝑢₀⁽⁰⁾𝑢₁⁽²⁾𝑥₁,𝑥₂−2𝑢₀⁽¹⁾𝑥₁𝑢₁⁽¹⁾𝑥₂,𝐿𝑢₂⁽³⁾=−2𝑢₀⁽¹⁾𝑥₁𝑢₁⁽²⁾𝑥₂,𝑥₃. (13) The zero boundary condition is applied to all kernels. The initial condition in (1) will be applied only to 𝑢₀⁽⁰⁾ while all other kernels will have zero initial conditions. This will result in some PDEs in (13) with zero initial and boundary conditions in addition to zero forcing terms. These kernels will be zeros. So, we will have 𝑢₀⁽²⁾=𝑢₀⁽³⁾=𝑢₁⁽³⁾=0. (14) The expectation and variance for the 𝑚th order 𝑁𝐶th correction solution will be obtained as follows [10]: 𝐸[𝑢(𝑡,𝑥)]=^𝑁𝐶_𝑖=0𝜆^𝑖𝑢_𝑖⁽⁰⁾,Var[𝑢(𝑡,𝑥)]=^𝑚_𝑗=1𝑗!_𝑅𝑗^𝑁𝐶_𝑖=0𝜆^𝑖𝑢_𝑖^(𝑗)2𝑑𝜏_𝑗. (15)(16) The same numerical solution described above (FVM + Crank-Nicolson) for WHE can be used to solve the deterministic system in (13) but, in this case, no iterations are required. The equations in (13) are decoupled and the right-hand side for each kernel 𝑢_𝑖^(𝑗) will be in terms of already (previously) computed kernels. For example, when solving for the kernel 𝑢₂⁽⁰⁾, the right-hand side is in terms of the kernels 𝑢₀⁽⁰⁾, 𝑢₁⁽⁰⁾, and 𝑢₁⁽¹⁾ which are already computed before.

5. The Analytical Solution

The deterministic system in (13) resulting from the WHEP technique can be solved analytically. In the current work Mathematica 8.0 will be used to solve analytically the WHEP system (13) as described below.

Consider the following deterministic model equation for any 𝑣=𝑢_𝑖^(𝑗);0≤𝑗≤𝑚,0≤𝑖≤𝑁𝐶: 𝐿(𝑣((𝑥);()∈()×(),𝑣()()𝑣()(𝑥).𝑡,𝑥))=𝑔𝑡,𝑥0,∞0,𝑡,0=𝑣𝑡,=0,0,𝑥=𝜑 (17) Apply the eigenfunction expansion technique [21] to get the following solution [11]: 𝑣(𝑡,𝑥)=^∞_𝑘=0𝑇_𝑘𝑒^{−(𝑘𝜋/)2𝑡}sin𝑘𝜋𝑥+^∞_𝑘=0𝐼_𝑘(𝑡)sin𝑘𝜋𝑥, (18) where 𝑇_𝑘=2₀𝜑(𝑥)sin𝑘𝜋𝐼𝑥𝑑𝑥,_𝑘(𝑡)=^𝑡₀𝑒^{−(𝑘𝜋/)2(𝑡−𝜏)}𝐹_𝑘(𝜏)𝑑𝜏, (19) where 𝐹_𝑘(𝑡)=2₀𝑔(𝑥)sin𝑘𝜋𝑥𝑑𝑥. (20) The two functions 𝜑(𝑥) and 𝑔(𝑥) for each kernel 𝑢_𝑖^(𝑗) areas are given in (13). This model solution can be written directly into Mathematica for each kernel. The expectation and variance of the solution are obtained using the same formulae in (16).

The advantage in using the analytical solution is that there are no restrictions on the solution convergence; that is, the solution can be always obtained and there are no limitations even on the values of the nonlinearity strength 𝜆.

6. Monte-Carlo Simulations (MCS)

Monte-Carlo simulations can be done by sampling the white-noise term 𝑁(𝑥;𝑤) and use the generated samples to solve the heat equation for each sample, that is, solving 𝐿[𝑢(𝑡,𝑥)]=−𝜆𝑢²+𝑓(𝑡)𝑁_𝑖(𝑥), (21) where 𝑁_𝑖(𝑥) is the 𝑖th sample of 𝑁(𝑥;𝑤). In the current work, the same numerical technique described above (FVM + Crank-Nicolson scheme) will be used to solve (21). The nonlinear term in the right-hand side of (21) will be computed explicitly. The expectation and variance are then obtained from the resulting set of solutions.

An arbitrary sample 𝑁_𝑖(𝑥) of the white noise 𝑁(𝑥;𝑤) can be generated with the well-known spectral representation [12]: 𝑁_𝑖(𝑥)=^𝑀𝜔_𝑗=14𝑆₀Δ𝜔cos𝜔_𝑗𝑥+𝜑_𝑗;1≤𝑖≤𝑛MC, (22) where 𝑛MC is the number of Monte-Carlo samples, Δ𝜔 is a constant step on the frequency axis, 𝜔_𝑗 are 𝑀_𝜔 equally spaced frequencies, and 𝜑_𝑗 are 𝑀_𝜔 random phases uniformly distributed in the interval [0,2𝜋]. In the current work, we will use Δ𝜔=0.05 and 𝑀_𝜔=500 to generate the samples.

7. Results and Comparisons

The initial solution is taken as 𝑢(0,𝑥)=𝑔(𝑥)=𝑥(−𝑥) and the deterministic envelope function 𝑓(𝑡)will be assumed unity.

The four techniques described above are tested for different values of the nonlinearity strength parameter 𝜆. In our case and with the other variables fixed, the parameter 𝜆 controls the convergence of the numerical algorithms (WHE, WHEP, and MCS). This convergence condition represents also the contraction property of the operator 𝐿⁻¹ in (9).

In the current work, it was found numerically that the maximum value of 𝜆 that results in convergent solution is around 10.5. Higher values of 𝜆 will cause divergence. The analytical solution does not suffer from this problem. It can be used for any value of 𝜆 but on the other hand it is not easy to be computed and requires more time and memory storage and sometimes may not be reachable at all specially for higher-order solutions.

Figure 1 shows the steady-state solution along the 𝑥-axis for 𝜆= 0.5 using the four algorithms. The numerical WHEP and numerical WHE give approximately the same solution. The MCS requires around 10⁶ samples to give convergent solution. The analytic solution gives a decayed solution compared with other solutions. As it is shown in the figure, the steady-state mean solution converges to a small negative value and not to zero due to the continuous excitation noise term.

(a)

(b)

Figure 1 

Steady-state mean (a) and variance (b) along 𝑥-axis for 𝜆 = 0.5 using the four different techniques.

A 3D plot for the expectation and variance is given in Figure 2 for 𝜆= 0.5. As it is shown in the figure, the mean solution decays with the time until it reaches a steady-state value that depends on the nonlinearity strength 𝜆. The variance starts from zero and increases with the time until it reaches its steady-state value.

(a)

(b)

Figure 2 

Solution mean (a) and variance (b) for 𝜆= 0.5.

Figure 3 shows the evolution of the mean solution with the time until steady state is reached at 𝑥=0.5 and 𝜆=0.5. The WHE mean solution converges to −0.00069 while the mean obtained from 1000,000 Monte-Carlo simulations converges to −0.0007 as shown in the figure. For the same parameters, the steady-state variance from WHE converges to 0.0189 while it converges to 0.0187 with Monte-Carlo simulations.

(a)

(b)

Figure 3 

Variation of the mean solution (a) and the variance (b) with the time using WHE and 106 MCS simulations at x = 0.5 and 𝜆= 0.5.

The steady-state values of the mean and variance vary with the number of Monte-Carlo samples. Figure 4 shows the variation of the mean and variance steady states with different number of Monte-Carlo simulations. As it is shown in the figure, the steady-state variance converges rapidly to approximately 0.018 with the number of Monte-Carlo samples while the steady-state mean value oscillates around −0.0007 with the number of Monte-Carlo samples.

(a)

(b)

Figure 4 

Steady-state mean (a) and steady-state variance (b) with the number of Monte-Carlo samples at x = 0.5 and 𝜆= 0.5.

Figure 5 shows the mean solution for different values of the nonlinearity strength 𝜆. As we can notice from the figure, increase in the nonlinearity strength 𝜆 will cause the steady-state value of the mean solution to increase in the negative direction.

Figure 5 

Mean solution for different values of the nonlinearity strength 𝜆.

The convergence rate depends on the value of the nonlinearity strength 𝜆. Figure 6 shows the convergence rates for different values of 𝜆 using WHE technique. As it is shown in the figure, the convergence rate slows down as 𝜆 increases. As the value of 𝜆 approaches 10, the convergence rate becomes very slow and for 𝜆>10.5, the solution starts to diverge when using Pickard’s iteration.

Figure 6 

Convergence rate for different values of the nonlinearity strength 𝜆.

The advantage of higher-order solutions is to estimate the non-Gaussian part of the solution especially in case of existence of nonlinearities. Figure 7 shows comparison between the total variance (Gaussian and non-Gaussian) and the non-Gaussian variance of the solution at x = 0.5 for 𝜆= 0.5. Figure 8 shows the percentage of the non-Gaussian variance to the total variance as the nonlinearity strength 𝜆 increases. For 𝜆= 10.0, the non-Gaussian variance is around 15% of the total variance. This means that, for higher nonlinearity strengths, the non-Gaussian effect cannot be neglected and hence higher-order WHE solutions are necessary.

(a)

(b)

Figure 7 

Total variance (a) and the non-Gaussian variance (b) for different values of the nonlinearity strength 𝜆.

Figure 8 

The percentage of the non-Gaussian variance to the total variance with the nonlinearity strength 𝜆.

The numerical WHEP algorithm consumes 33 seconds to get the solution while the numerical WHE requires 83 seconds to converge with residual of 10⁻⁶. The MCS consumes 11.3 seconds per 1000 samples in our case. This emphasizes the fact that MCS is not an efficient algorithm and it can be used only for validation of other algorithms. The analytical solution is not efficient as well and it consumes around 900 seconds.

Numerical tests show that numerical WHE and WHEP algorithms are more efficient compared with the analytical solution and, of course, with the time-consuming MCS. Numerical WHEP technique is faster than numerical WHE algorithm because we do not have iterations in the WHEP algorithm. On the other hand, WHEP algorithm is memory consuming as there is large number of kernels to be computed. The WHEP has (𝑚+1)(𝑁𝐶+1) kernels while WHE has only (𝑚+1) kernels. For higher orders and higher corrections the WHEP algorithm will not be practical as it may not fit in the computer memory. We can conclude from these discussions that the numerical WHE algorithm is the most practical algorithm among other algorithms to solve the SPDEs with higher-order approximations.

8. Conclusions

In the current work, the stochastic heat equation with nonlinearities is considered. The stochastic equation is converted to an equivalent deterministic system by applying the WHE technique. The resulting deterministic coupled system is solved numerically using FVM and the Crack-Nicolson scheme. The convergence condition of the solution was discussed and estimated numerically. The perturbation technique is applied to the coupled deterministic equivalent system to decouple it. The same numerical scheme is used to solve the decoupled system. The solution is validated by comparing with the analytical solution and the MCS of the original stochastic equation. The comparisons show the efficiency in using WHE in solving the stochastic PDEs. The results also show the importance of higher-order WHE solution especially in case of nonlinearities. The Gaussian and non-Gaussian solutions are compared and it was concluded that the non-Gaussian effects cannot be neglected and hence there is a need for higher-order solutions.

Conflict of Interests

The authors declare that they have no conflict of interests regarding to the publication of this paper.