Journal of Applied Mathematics

Volume 2014 (2014), Article ID 843714, 9 pages

http://dx.doi.org/10.1155/2014/843714

Research Article

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

Mohamed A. El-Beltagy1,2 and Noha A. Al-Mulla³

¹Electrical & Computer Engineering Department, Engineering Faculty, Effat University, Jeddah 21478, Saudi Arabia

²Engineering Mathematics & Physics Department, Engineering Faculty, Cairo University, Giza 12221, Egypt

³Mathematics Department, Science College, University of Dammam, Dammam 31451, Saudi Arabia

Received 5 April 2014; Accepted 8 June 2014; Published 24 June 2014

Academic Editor: Song Cen

Copyright © 2014 Mohamed A. El-Beltagy and Noha A. Al-Mulla. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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   in [17], (1) has a unique finite mean solution if the functions and are globally Lipschitz and if the initial condition . 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

In the current work and without loss of generality, a unit thermal diffusivity () and second degree nonlinear losses are assumed. Now, (1) can be rewritten as

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]

where is the deterministic kernel of , , and is a -dimensional integral over the disposable variables . The functional is the 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 kernels ; that is, (2)substitute the truncated expansion into  (3),(3)use the multinomial theorem to expand the nonlinear term ,(4)multiply by and then apply the ensemble average.

This will result in () equations in the deterministic kernels as

The expectations are computed as . 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 equals one when 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 In case of , 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 , the expectation and variance of the solution are obtained as [10]

In case of quadratic nonlinearity () and 3rd order approximation (), (5) will expand to the following deterministic system:

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 :

Pickard’s iteration takes the form

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:

where is the kernel at time level and space subdivision . The factor is defined as , and is defined as

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

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 , apply the perturbation technique up to corrections; that is, substitute in (7).(6)Compare the coefficients of in both sides to get equations for each kernel .

This will result in 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 will be

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

The expectation and variance for the order correction solution will be obtained as follows [10]:

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 :

Apply the eigenfunction expansion technique [21] to get the following solution [11]:

where

where

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

where is the 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]:

where 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 . In the current work, we will use and to generate the samples.

7. Results and Comparisons

The initial solution is taken as 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 and . 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 10⁶ 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 , 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 kernels while WHE has only 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.

References

J. Cannon, The One-Dimensional Heat Equation, Cambridge University Press, Cambridge, UK, 1984.

G. A. Bécus, “Random generalized solutions to the heat equation,” Journal of Mathematical Analysis and Applications, vol. 60, no. 1, pp. 93–102, 1977.

R. Marcus, “Parabolic Ito equations with monotone nonlinearities,” Journal of Functional Analysis, vol. 29, no. 3, pp. 275–286, 1978.

R. Manthey, “Weak convergence of solutions of the heat equation with Gaussian noise,” Mathematische Nachrichten, vol. 123, pp. 157–168, 1985.

R. Manthey, “Existence and uniqueness of a solution of a reaction-diffusion equation with polynomial nonlinearity and white noise disturbance,” Mathematische Nachrichten, vol. 125, pp. 121–133, 1986.

G. Jetschke and R. Manthey, “On a linear reaction diffusion equation with white noise on the whole real line,” Forschungsergebnisse N/86/17, Friedrich Schiller University of Jena, Jena, Germany, 1986.

M. Foondun and D. Khoshnevisan, “On the stochastic heat equation with spatially-colored random forcing,” Transactions of the American Mathematical Society, vol. 365, no. 1, pp. 409–458, 2013.

Z. Brzeźniak, E. Hausenblas, and P. Razafimandimby, “Martingale solutions for Stochastic Equation of Reaction Diffusion Type driven by Lévy noise or Poisson random measure,” http://arxiv.org/abs/1010.5933v2

E. F. Abdel Gawad and M. A. El Tawil, “General stochastic oscillatory systems,” Applied Mathematical Modelling, vol. 17, no. 6, pp. 329–335, 1993.

M. A. El-Beltagy and M. A. El-Tawil, “Toward a solution of a class of non-linear stochastic perturbed PDEs using automated WHEP algorithm,” Applied Mathematical Modelling: Simulation and Computation for Engineering and Environmental Systems, vol. 37, no. 12-13, pp. 7174–7192, 2013.

M. A. El-Tawil and N. A. Al-Mulla, “The solution of non-linear diffusion equation under stochastic nonhomogeneity using symbolic {WHEP} and Pickard algorithms,” in Transactions on Computational Science VII, vol. 5890 of Lecture Notes in Computer Science, pp. 75–100, Springer, Berlin, Germany, 2010.

M. El-Beltagy and A. Al-Johani, “Stochastic response of duffing oscillator with fractional or variable-order damping,” Journal of Fractional Calculus and Applications, vol. 4, no. 2, pp. 357–366, 2013.

A. Al-Johani and M. El-Beltagy, “Numerical solution of stochastic nonlinear differential equations using wiener-hermite expansion,” in Proceedings of the 11th International Conference of Numerical Analysis and Applied Mathematics (ICNAAM '13), vol. 1558 of AIP Conference Proceedings, pp. 2099–2102, 2013.

M. Hamed, M. El-Tawil, B. El-Desouky, and M. El-Beltagy, “Solution of nonlinear stochastic Langevin's equation using WHEP, Pickard and HPM methods,” Applied Mathematics, vol. 5, no. 3, pp. 398–412, 2014.

A. Fareed, H. El-Zoheiry, M. El-Tawil, M. El-Beltagy, and H. Hassan, “Solving nonlinear stochastic diffusion models with nonlinear losses using the homotopy analysis method,” Applied Mathematics, vol. 5, no. 1, pp. 115–127, 2014.

M. El-Beltagy and A. Al-Johani, “Numerical approximation of higher-order solutions of the quadratic nonlinear stochastic oscillatory equation using WHEP technique,” Journal of Applied Mathematics, vol. 2013, Article ID 685137, 21 pages, 2013.

R. Dalang, D. Khoshnevisan, C. Mueller, D. Nualart, and Y. Xiao, A Minicourse on Stochastic Partial Differential Equations, Springer, Salt Lake City, Utah, USA, 2006.

X. Bardina, M. Jolis, and L. Quer-Sardanyons, “Weak convergence for the stochastic heat equation driven by Gaussian white noise,” Electronic Journal of Probability, vol. 15, no. 39, pp. 1267–1295, 2010.

V. Berinde, “Iterative approximation of fixed points for pseudo-contractive operators,” Seminar on Fixed Point Theory, vol. 3, pp. 209–215, 2002.

M. Giles and R. Carter, “Convergence analysis of Crank-Nicolson and Rannacher time-marching,” Journal of Computational Finance, vol. 9, no. 4, pp. 89–112, 2006.

S. J. Farlow, Partial Differential Equations for Scientists and Engineers, John Wiley & Sons, New York, NY, USA, 1982.