Abstract
This paper deals with a research question raised by Jentzen and Röckner (A Milstein scheme for SPDEs, arXiv:1001.2751v4 (2012)), whether the exponential term in their introduced scheme can be replaced by a simpler mollifier. This replacement can lead to more simplification and computational reduction in simulation. So, in this paper, we essentially replace the exponential term with a Padé approximation of order 1 and denote the resulting scheme by simplified Milstein scheme. The convergence analysis for this scheme is carried out and it is shown that even with this replacement the order of convergence is maintained, while the resulting scheme is easier to implement and slightly more efficient computationally. Some numerical tests are given that confirm the order of accuracy and also computational cost reduction.
1. Introduction
Many models in engineering, physics, complex phenomena, and so forth are described by stochastic partial differential equations (SPDEs); for example, see [1–6]. Since the exact solutions of these equations are rarely known, the numerical analysis of SPDEs has been recently the subject of many papers; for example, see [7–14], for more detailed discussion on this topic and many examples in applied sciences. In this paper, we consider strong approximation (see [4, Section 9.3]) of SPDEs of evolutionary type. To demonstrate the results of this paper clearly, we focus on the following example of SPDE: with initial condition and Dirichlet boundary conditions for all and , where . Let be a probability space with a normal filtration , let be the -Hilbert space of equivalence classes of Lebesgue square integrable functions from to , and let be two appropriate smooth and regular functions with globally bounded derivatives. Let be a standard -Wiener process with regard to , with a trace class operator and with being a smooth function. The covariance operator has orthonormal basis , , of eigenfunctions with summable eigenvalues . Under the previous assumption, the SPDE (1) has a unique mild solution. Specifically, there exists an up to unpredictable unique adapted stochastic process with continuous sample path which satisfies for all , where is the Laplacian with Dirichlet boundary conditions times the constant and where and are given by and for all , and all where with for all is the image -Hilbert space of (see [15, Appendix C]); note that and commutate in our example SPDE (2). Now we are concerned about the strong approximation of the SPDE (1). More formally we want to compute numerical approximation such that holds for a given precision with the least possible computational effort. To simulate the numerical approximation on a computer, one has to discretize both the time interval and the infinite dimensional space . In this paper we consider spectral Galerkin for spatial discretization and difference method for temporal discretization. A simple full discretization for (1) is the linear implicit Euler method combined with spectral Galerkin method which is given by , and all , with and , where is a bounded linear operator such that with for all , , and , and the finite dimensional Wiener processes , , are given by for all , , and . According to the analysis of [16], for method (4) with , there exist real numbers that for small holds for all . This means that the method has overall convergence (for a real number , we write for the convergence order if the convergence order is higher in order than for all arbitrary small ). In [16] Jentzen and Rockner proposed an infinite dimensional analog of Milstein type scheme for (1) given by and for all and . Here we use the notations , , and for all and all functions , . Method (7) gives a break of complexity of the numerical approximation of nonlinear SPDE with multiplicative trace class noise. More precisely, it is shown in [16] that time steps in contrast to time steps for the linear implicit Euler scheme (4) are required to achieve (6). That is the Milstein type scheme (7) with time steps guarantees that for real numbers , , such that holds for all . Thus the scheme has the overall convergence order of . Consequently scheme (7) increases the overall convergence order from to . As mentioned before, in this paper essentially the exponential term in the Milstein type scheme [16] is replaced by a first order approximation which makes the scheme easier to implement and slightly more efficient computationally while preserving the order of convergence. The analysis and implementation will be carried out as follows. In Section 2 the required setting and assumptions are formulated. In Section 3 the simplified Milstein scheme is formulated. In Section 4 we state and prove the main result of this section concerning the convergence of the simplified Milstein scheme. Finally in Section 5 numerical example for a stochastic reaction diffusion equation is presented to show numerically the order of convergence and computational costs. The numerical simulations will be carried out in MATLAB environment on a PC with CPU 2.66 GHz.
2. Setting and Assumptions
Throughout this paper suppose that the setting and following assumptions are fulfilled. Fix . Let be a probability space with a normal filtration and let and be two separable -Hilbert spaces. Moreover, let be a standard -Wiener process with respect to , with a trace class operator .
Assumption 1 (linear operator ). Let be a linear operator such that
for every with .
Here is a family of real numbers with and is an orthonormal basis of . By equipped with the norm for all , , we denote the -Hilbert space of domains of fractional powers of the linear operator .
Assumption 2 (drift term ). Let be a real number and let be a globally Lipschitz continuous; that is, and , ; in addition
Assumption 3 (diffusion term ). Let be a globally Lipschitz continuous mapping and twice continuously Frechet differentiable mapping with and . In addition , with , , and is a real number such that and
hold for all and . Additionally, let the bilinear Hilbert-Schmidt operator be symmetric for all . Note that the operator , given by
for all , is a bilinear Hilbert-Schmidt operator in for all .
The assumed symmetry of thus reads as [16, Remark 1].
Assumption 4 (initial value ). Let be an -measurable mapping with .
Proposition 5 (existence of the mild solution). Let . Then under Assumptions 1–4, there exists an up to modifications unique predictable stochastic process which fulfills , for all ; moreover, we have
Proposition 5 immediately follows from Theorem 1 in [17].
3. The Proposed Simplified Milstein Scheme
We construct the simplified Milstein scheme for nonlinear stochastic partial differential equations. For this work first we use Taylor formula in Banach space for coefficients and for the problem (2). More formally using and for small shows for small . Using the approximation for small gives We then substitute for small to obtain where . Combining the temporal approximation (19) and spatial discretization in (4) suggests the numerical scheme given by and for all and all , where . The difficulty in this formula is working with the term corresponding to the double integral. As suggested by Jentzen and Rockner (see [16, Subsection 6.7]), this double integral can be replaced by By using (21), the numerical scheme (20) thus reduces to where and for all , .
For the simplified Milstein scheme (22) applied to (1), the main result of this paper, that is, Theorem 7, will show that with Similar to scheme (7), the numerical method (22) can be simulated quite easily. The term in (22) can be computed once in advance for which computational operations are needed. With the term at hand, further computational operations and independent standard normal random variables are needed to compute from by using fast Fourier transform. Therefore, since time steps are used, computational operations and random variables are required to obtain . The logarithmic term in arises from fast Fourier transform computations, due to the nonlinearities of and . Taking into account the convergence order in (23), one can show that scheme (22) shares the same overall convergence order of , which is greater than the overall convergence order of the linear implicit scheme (4). We then take a more closer look at schemes (7) and (22) at each step. It is obvious that the Milstein scheme (7) requires evaluation of exponential term, while the simplified Milstein scheme needs one simple mollifier instead of exponential term. The CPU time for one path simulation by the simplified Milstein scheme (22) applied to (1) is less than that for (7). For example, see Table 1; for , one path simulation of the simplified Milstein scheme (22) requires 834.842365 CPU seconds, while this simulation by Milstein scheme (7) needs 842.699400 CPU seconds. This difference is due to the fact that evaluation of the exponential term takes more time than that of the simple mollifier term. A natural question thus arises on whether such substitution can maintain the high convergence order of (7). In this paper we investigate this issue and prove that the simplified Milstein scheme maintains the expected order of convergence.
4. Convergence Analysis
Let be an orthonormal basis consisting of the eigenfunctions of , and let be their corresponding eigenvalues with as a trace class operator; that is, for all . We define the linear projection operator Furthermore, we define Wiener processes by for all , , and . Let , , and be the time discretization step, and let be the -Hilbert space of equivalence classes of -measurable and Lebesgue square integrable functions from to with the scalar product and the norm for . Now we start our investigation to analyze the proposed method for SPDE fulfilling Assumptions 1–4. Based on (1) we then consider and , , as the orthonormal basis of , which satisfy For the drift term to fulfill Assumption 2, let be a continuously differentiable function with and Then, the operator given by , for and , satisfies Assumption 2. For the diffusion term to satisfy Assumption 3, we consider to be a twice continuously differentiable function with and also for all and some given . Let be the operator for all .
It has been shown in [16] that fulfills Assumption 3. For the initial value to satisfy Assumption 4, we assume that is a twice continuously differentiable function with . Then the mapping given by for all fulfills Assumption 4 for all . With the above setting, the SPDE (1) reduces to with and for , , and .
Moreover, we define a family for all , and all , and we consider the mappings , by .
Using these notations, the SPDE (30) can be rewritten as with and for and .
Scheme (21)-(22) applied to the SPDE (30) reduces to where and , , . Therefore the numerical method (32) satisfies where and for all , . The following inequalities are classical and one can easily prove them by using the spectral decomposition of [1]: To give the order of the convergence for the simplified Milstein approximation of the evolution equation, we need the following version of the Gronwall lemma.
Lemma 6. Let be two sequences of nonnegative numbers such that and such that there exists a positive constant such that then
Proof. By Mathematical induction with respect to using .
From the above lemma we can deduce that The main result of this section is stated below.
Theorem 7. Let , , and . Suppose that is the solution of (2) on . Let Assumptions 1–4 hold, and let be the numerical approximations obtained by scheme (33). Then there exists a positive constant such that
Proof. To start the proof, we first note that the exact solution of SPDE (2) satisfies
where . In particular, (43) shows
Let
where
For the spatial discretization error , we have
the real numbers are given by (see [16])
For the with respect to (33) and (44), we have
where and and are, respectively, the other terms under operator. From (36), the first term of (49) can be easily estimated by
Let denote a constant which may depend on , or . We now treat the second term of (49)
According to Assumption 2 and the fact that for all , we get
For the first term on the right-hand side of inequality (52), we have
To estimate the second term of (52), one should note that
and from (35) and (37), we have
Therefore, from (35), (38), and (36), we have
thus we have
Therefore, by taking the expectation of (52) to the power of 2 and using (53) and (57), we get
in which the summation can be estimated as
For the second term of (51), Proposition 5 and Assumption 2 lead to
For the third term of (51), by using Assumption 2, we have
which leads to
where we have used the Minkowski inequality in (62).
Therefore, from (58), (60), and (62), we get
Now for the last term of (49), we obtain
where
Using the fact that for all , for , we have
and therefore
Thus
which means
Therefore, from (13), we get
which leads to
For , by Proposition 5, it is seen that
in which
Now for the first term on the right-hand side of (73), we need to estimate
and then by (34) and (36), we obtain
Similarly from (34), (35), we obtain
where is such that
which is possible since
Therefore,
For the second term on the right-hand side of (73), we can write
from which we get
Thus from (73), after replacing (79) and (81) in (72), we obtain
which gives
For the third term , we have
This implies that
Finally for , the last term of (65), we should recall
Using the fact that
for all , , , and and using the inequality
for all , , for , we obtain
To estimate the first term of (89), we first approximate
for all , with and all . More precisely, with respect to (88), we have
By using (71), we have
This implies that
for all , with , and all . Therefore, we obtain
and hence
for all , with .
Therefore,
and hence from (95), we get
And for the second term of (89), we have
Therefore, from (71), (83), (85), (89), (97), and (98), we obtain
Hence from (50), (63), and (100), we obtain
Now we take an integer and use the Holder inequality for the two summations in the last estimation to get
Therefore, with using (97) and (98), we obtain
We then have
Hence, we conclude from (41) that
Finally, with respect to (47) and (105), we obtain
which completes the proof of the theorem.
5. Simulation Results
In this section we consider SPDE (1) and solve it by numerical scheme (22). More formally, let and be given by for all and suppose that are given by and for all , . The SPDE (1) reduces to with and for and . We also assume that the SPDE (107) should be solved with a precision of, say, two decimals, that is, with the precision in (3). To confirm numerically our theoretical founding in Theorem 7, we recall that for SPDE (107) there should exist some real number such that holds for each small . The overall convergence order of the linear implicit Euler method (4) is (see [16]), while the overall convergence of the simplified Milstein scheme (22) and Milstein scheme (7) is . In Figure 1 the approximation error in the sense of (6) of the linear implicit Euler approximation , obtained by (4), of the approximation , obtained by Milstein scheme (7), and of the approximation , obtained by simplified Milstein scheme (22), is plotted against the precise number of independent standard normal random variables that is needed to compute the corresponding approximation for on a log-log scale. Figure 1 confirms the order of convergence of our scheme and compares with the other two schemes. Besides, the simplified Milstein scheme (22) and the Milstein scheme (7) produce nearly the same approximation errors. Numerical results also show that the linear simplified Milstein scheme (21) and the Milstein type scheme (7) are much computationally effective than the linear implicit Euler scheme (4). To simulate one path , one needs to generate independent normal random variables, but this amount for simulation of and reduces to . From the numerical results reported in Table 1 and Figure 1, we conclude that the simplified Milstein scheme is more effective than implicit Euler method and slightly better than Milstein scheme.
6. Conclusions
A simplified Milstein scheme for solving stochastic partial differential equations of the form (1) with multiplicative trace class noise was theoretically and numerically investigated. This scheme has advantages to some other methods such as linear implicit Euler and Milstein schemes. We have shown the convergence of this method under the stated conditions and then we have illustrated the effectiveness of the simplified Milstein scheme numerically.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.