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

Authors: Mohamed A. El-Beltagy; Noha A. Al-Mulla
Journal: Journal of Applied Mathematics
 (2014)

---

## 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.


---

## Body

## 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 Section2. 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(1)∂u∂t-α∂2u∂x2=-λun+f(t)N(x;w);1111111111(t,x)∈(0,∞)×(0,l),u(t,0)=u(t,l)=0,u(0,x)=φ(x),
where α is the thermal diffusivity and ω is the random outcome of a triple probability space (Ω,χ,P) in which Ω is a sample space, χ is σ-field associated with Ω, and P is a probability measure. The nonlinear losses term un is scaled (strengthened) by a deterministic scalar λ. The term N(x;w) is a space white-noise excitation term that has the property E[N(x;w),N(x1;w)]=δ(x-x1). The function f(t) is a deterministic envelope function. The initial condition φ(x) will be assumed to be deterministic. The nonlinear losses term un may account for the radiative heat loss which depends on the excess heat u compared with the surroundings. At low excess temperatures, the radiative loss may be approximated linearly (proportional to u) while at high excess temperature, the Stefan-Boltzmann law gives a net radiative heat loss proportional to u4.According to theorem4.1 in [17], (1) has a unique finite mean solution if the functions λun and f(t) are globally Lipschitz and if the initial condition φ(x)∈C([0,l]). 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 operatorL as
(2)L=∂∂t-α∂2∂x2.
In the current work and without loss of generality, a unit thermal diffusivity (α=1) and second degree nonlinear losses (n=2) are assumed. Now, (1) can be rewritten as
(3)L[u(t,x;w)]=-λu2+f(t)N(x;w),
with the initial and boundary conditions described in (1).


## 3. Description of the WHE Technique

In the WHE technique, the stochastic response functionu(t,x;w) is expanded as [10]
(4)u(t,x;w)=∑k=0∞∫Rku(k)(t,x;x1,x2,…,xk)111111111111111×H(k)(x1,x2,…,xk)dτk,
where u(k)(t,x;x1,x2,…,xk) is the kth deterministic kernel of u(t,x;w), dτk=dx1dx2⋯dxk, and ∫Rk‍ is a k-dimensional integral over the disposable variables x1,x2,…,xk. The functional H(k)(x1,x2,…,xk) is the kth order Wiener-Hermite functional. The WHE technique for general nonlinear exponent (n) and general order (m) can be summarized in the following steps [10]. For the sake of brevity, the space and time variables (t,x) and ω will be omitted and only disposable variables x1,x2,x3,… are written: (1)
truncate the expansion (4) to contain only m+1 kernels u(j);0≤j≤m; that is, u(t,x;w)=u(0)(t)+∑j=1m∫Rju(j)H(j)dτj,(2)
substitute the truncated expansion into  (3),(3)
use the multinomial theorem to expand the nonlinear termun;n=2,(4)
multiply byH(j);0≤j≤m and then apply the ensemble average.This will result in (m+1) equations in the deterministic kernels u(j);0≤j≤m as
(5)L(u(j))=-λj!∑gcg∫Rz(∏i=0m[u(i)]kgi)Egjdτz+δj1j!f(t)δ(x-x1);111111110≤j≤m.
The expectations Egj are computed as Egj=〈H(j)∏i=0m(H(i))kgj〉. It was explained in [10] how to get Egj in terms of the Dirac delta functions and then use them to reduce the order of integration. The Kronecker delta function δj1 equals one when j=1 and zero otherwise. The counter g in the summation in the right-hand side of (5) runs over all the (n+mn) combinations of the positive integers kg0,kg1,…,kgm such that ∑i=0mkgi=n. In case of m=1, the Gaussian part of the solution is obtained. Larger values of m will account for the non-Gaussian part of the solution.After solving (5) for the kernels u(j);0≤j≤m, the expectation and variance of the solution are obtained as [10]
(6)E[u(t,x;w)]=u(0),Var⁡[u(t,x;w)]=∑k=1m(k!)∫Rk[u(k)]2dτk.
In case of quadratic nonlinearity (n=2) and 3rd order approximation (m=3), (5) will expand to the following deterministic system:
(7)L(u(0))=-λ[u(0)]2-λ∫R[u(1)(x1)]2dx111-2λ∫R2[u(2)(x1,x2)]2dx1dx211-6λ∫R3[u(3)(x1,x2,x3)]2dx1dx2dx3,L(u(1))=f(t)δ(x-x1)-2λu(0)u(1)(x1)11-4λu(1)(x2)u(2)(x1,x2)dx211-12λ∫R2u(2)(x2,x3)u(3)(x1,x2,x3)dx2dx3,L(u(2))=-2λu(0)u(2)(x1,x2)-λu(1)(x1)u(1)(x2)11-6λ∫Ru(1)(x3)u(3)(x1,x2,x3)dx311-4λ∫Ru(2)(x1,x3)u(2)(x2,x3)dx311-18λ∫R2u(3)(x1,x3,x4)u(3)(x2,x3,x4)dx3dx4,L(u(3))=-2λu(0)u(3)(x1,x2,x3)-2λu(1)(x1)u(2)(x2,x3)11-4λ∫Ru(2)(x1,x4)u(3)(x2,x3,x4)dx411-4λ∫Ru(2)(x2,x4)u(3)(x1,x3,x4)dx411-4λ∫Ru(2)(x3,x4)u(3)(x1,x2,x4)dx4.
The zero boundary condition is applied to all kernels in (7). The initial condition in (1) will be applied only to u(0) 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 v=u(j);0≤j≤m:
(8)L(v)=F(v,u(0),…,u(m));v(0,x;x1,…,xj)=v0.
Pickard’s iteration takes the form
(9)vk+1=L-1[F(v,u(0),…,u(m))|k],
where k 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 L-1 in (9). If the operator L-1 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 L-1 is a strict contraction or satisfies a certain generalized strong contraction condition. In our case, the contraction condition of L-1 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 Δt and Δx 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:
(10)-0.5rvi-1k+1+(1+r)vik+1-0.5rvi+1k+1=0.5rvi-1k+(1-r)vik-0.5rvi+1k+ΔtΔxF-ik,
where vik is the kernel v at time level k and space subdivision i. The factor r is defined as r=Δt/(Δx)2, and F-ik is defined as
(11)F-ik=∫xi-Δx/2xi+Δx/2Fk(v,u(0),…,u(m))dx.
The above described scheme in (10) has a truncation error of O(Δt2,Δx2) [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
(12)Residual=∑j=0m∥uk+1(j)-uk(j)∥.


## 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 NC. In the WHEP technique, the following two steps will be added to the above 4 steps of the WHE technique.(5)
For each kernelu(j);0≤j≤m, apply the perturbation technique up to NC corrections; that is, substitute u(j)=∑i=0NCλiui(j) in (7).(6)
Compare the coefficients ofλk;0≤k≤NC in both sides to get NC+1 equations for each kernel u(j);0≤j≤m.This will result in (m+1)(NC+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 NC=2 will be
(13)L(u0(0))=0,L(u0(1))=f(t)δ(x-x1),L(u0(2))=0,L(u0(3))=0,L(u1(0))=-[u0(0)]2-∫R[u0(1)(x1)]2dx1,L(u1(1))=-2u0(0)u0(1)(x1),L(u1(2))=-u0(1)(x1)u0(1)(x2),L(u1(3))=0,L(u2(0))=-2u0(0)u1(0)-2∫Ru0(1)(x1)u1(1)(x1)dx1,L(u2(1))=-2u0(0)u1(1)(x1)-2u1(0)u0(1)(x1)-4∫Ru0(1)(x2)u1(2)(x1,x2)dx2,L(u2(2))=-2u0(0)u1(2)(x1,x2)-2u0(1)(x1)u1(1)(x2),L(u2(3))=-2u0(1)(x1)u1(2)(x2,x3).
The zero boundary condition is applied to all kernels. The initial condition in (1) will be applied only to u0(0) 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
(14)u0(2)=u0(3)=u1(3)=0.
The expectation and variance for the mth order NCth correction solution will be obtained as follows [10]:
(15)E[u(t,x)]=∑i=0NCλiui(0),(16)Var⁡[u(t,x)]=∑j=1m(j!)∫Rj(∑i=0NCλiui(j))2dτj.
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 ui(j) will be in terms of already (previously) computed kernels. For example, when solving for the kernel u2(0), the right-hand side is in terms of the kernels u0(0), u1(0), and u1(1) 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 anyv=ui(j);0≤j≤m,0≤i≤NC:
(17)L(v(t,x))=g(x);(t,x)∈(0,∞)×(0,l),v(t,0)=v(t,l)=0,v(0,x)=φ(x).
Apply the eigenfunction expansion technique [21] to get the following solution [11]:
(18)v(t,x)=∑k=0∞Tke-(kπ/l)2tsin(kπl)x+∑k=0∞Ik(t)sin(kπl)x,
where
(19)Tk=2l∫0lφ(x)sin(kπl)xdx,Ik(t)=∫0te-(kπ/l)2(t-τ)Fk(τ)dτ,
where
(20)Fk(t)=2l∫0lg(x)sin(kπl)xdx.
The two functions φ(x) and g(x) for each kernel ui(j) 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 termN(x;w) and use the generated samples to solve the heat equation for each sample, that is, solving
(21)L[u(t,x)]=-λu2+f(t)Ni(x),
where Ni(x) is the ith sample of N(x;w). 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 sampleNi(x) of the white noise N(x;w) can be generated with the well-known spectral representation [12]:
(22)Ni(x)=∑j=1Mω4S0Δωcos⁡(ωjx+φj);1≤i≤nMC,
where nMC is the number of Monte-Carlo samples, Δω is a constant step on the frequency axis, ωj are Mω equally spaced frequencies, and φj are Mω random phases uniformly distributed in the interval [0,2π]. In the current work, we will use Δω=0.05 and Mω=500 to generate the samples.


## 7. Results and Comparisons

The initial solution is taken asu(0,x)=g(x)=x(l-x) and the deterministic envelope function f(t)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 L-1 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.Figure1 shows the steady-state solution along the x-axis for λ= 0.5 using the four algorithms. The numerical WHEP and numerical WHE give approximately the same solution. The MCS requires around 106 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.Steady-state mean (a) and variance (b) alongx-axis for λ = 0.5 using the four different techniques.
(a)
(b)A 3D plot for the expectation and variance is given in Figure2 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.Solution mean (a) and variance (b) forλ= 0.5.
(a)
(b)Figure3 shows the evolution of the mean solution with the time until steady state is reached at x=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.Variation of the mean solution (a) and the variance (b) with the time using WHE and 106 MCS simulations atx = 0.5 and λ= 0.5.
(a)
(b)The steady-state values of the mean and variance vary with the number of Monte-Carlo samples. Figure4 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.Steady-state mean (a) and steady-state variance (b) with the number of Monte-Carlo samples atx = 0.5 and λ= 0.5.
(a)
(b)Figure5 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. Figure7 shows comparison between the total variance (Gaussian and non-Gaussian) and the non-Gaussian variance of the solution atx = 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.Total variance (a) and the non-Gaussian variance (b) for different values of the nonlinearity strengthλ.
(a)
(b)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−6. 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(m+1)(NC+1) kernels while WHE has only (m+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.



---

*Source: 843714-2014-06-24.xml*