Numerical solutions of stochastic PDEs driven by arbitrary type of noise

Abstract

So far the theory and numerical practice of stochastic partial differential equations (SPDEs) have dealt almost exclusively with Gaussian noise or Lévy noise. Recently, Mikulevicius and Rozovskii (Stoch Partial Differ Equ Anal Comput 4:319–360, 2016) proposed a distribution-free Skorokhod–Malliavin calculus framework that is based on generalized stochastic polynomial chaos expansion, and is compatible with arbitrary driving noise. In this paper, we conduct systematic investigation on numerical results of these newly developed distribution-free SPDEs, exhibiting the efficiency of truncated polynomial chaos solutions in approximating moments and distributions. We obtain an estimate for the mean square truncation error in the linear case. The theoretical convergence rate, also verified by numerical experiments, is exponential with respect to polynomial order and cubic with respect to number of random variables included.

Appendix A: Interaction coefficients \(B(\alpha ,\beta ,p)\)

We still assume that $\{\xi _k\}_{k=1}^{\infty }$ are i.i.d. random variables. Then the interaction coefficient $B(\alpha ,\beta ,p)$ can be decomposed into

$$\begin{aligned} B(\alpha ,\beta ,p)= & {} \frac{{\mathbb {E}}[\Phi _{\alpha }\Phi _{\beta }\Phi _{p}]}{\alpha !}\nonumber \\= & {} \prod _{k=1}^{\infty }\frac{{\mathbb {E}}[\varphi _{\alpha _k}(\xi _k)\varphi _{\beta _k}(\xi _k)\varphi _{p_k}(\xi _k)]}{\alpha _k!}:=\prod _{k=1}^{\infty }b(\alpha _k,\beta _k,p_k). \end{aligned}$$

(A.1)

It suffices to compute b(i, j, l) for any $i,j,l\ge 0$. According to orthogonality,

$$\begin{aligned} \varphi _j(\xi )\varphi _l(\xi )=\sum _{i=0}^{\infty }\frac{{\mathbb {E}}[\varphi _i(\xi )\varphi _j(\xi )\varphi _l(\xi )]}{i!}\varphi _i(\xi )=\sum _{i=0}^{\infty }b(i,j,l)\varphi _i(\xi ). \end{aligned}$$

(A.2)

Hence b(i, j, l) is the ith expansion coefficient of $\varphi _j(\xi )\varphi _l(\xi )$ in terms of $\{\varphi _n(\xi )\}_{n=0}^{\infty }$. In particular, for the three types of noises and corresponding orthogonal polynomials considered in Section 4, there are explicit formulas for these expansion coefficients.

• For Gaussian noise and Hermite chaos. $\varphi _n(\xi )=He_n(\xi )$. Since

  $$\begin{aligned} He_j(x)He_l(x)=\sum _{r=0}^{\min \{j,l\}}\frac{j!l!}{(j-r)!(l-r)!r!}He_{j+l-2r}(x), \end{aligned}$$

  (A.3)

  we have

  $$\begin{aligned} b(i,j,l)=\left\{ \begin{array}{ll} \frac{j!l!}{(j-r)!(l-r)!r!} &{}\quad \text {if }i=j+l-2r\text { and }r\le \min \{i,j\}\\ 0 &{}\quad \text {otherwise} \end{array}\right. . \end{aligned}$$

  (A.4)

• For uniform noise and Legendre chaos, $\varphi _n(\xi )=\sqrt{(2n+1)n!}L_n(\xi /\sqrt{3})$. Define

  $$\begin{aligned} \lambda _n:=\frac{\Gamma (n+1/2)}{n!\Gamma (1/2)}=\frac{\prod _{m=0}^{n-1}(m+1/2)}{n!}. \end{aligned}$$

  Then the expansion of $L_j(x)L_l(x)$ is

  $$\begin{aligned} L_j(x)L_l(x)=\sum _{r=0}^{\min \{j,l\}}\frac{2(j+l-2r)+1}{2(j+l-r)+1}\frac{\lambda _r\lambda _{i-r}\lambda _{j-r}}{\lambda _{i+j-r}}L_{j+l-2r}(x). \end{aligned}$$

  (A.5)

  Thus

  $$\begin{aligned} b(i,j,l)=\left\{ \begin{array}{ll} \frac{\sqrt{(2i+1)(2j+1)(2l+1)}}{2(j+l-r)+1}\sqrt{\frac{j!l!}{i!}}\frac{\lambda _r\lambda _{i-r}\lambda _{j-r}}{\lambda _{i+j-r}} &{}\quad \text {if }i=j+l-2r\text { and }r\le \min \{i,j\}\\ 0 &{}\quad \text {otherwise} \end{array}\right. . \end{aligned}$$

  (A.6)

• For Beta$(\frac{1}{2},\frac{1}{2})$ noise and Chebyshev chaos, $\varphi _n(\xi )=\sqrt{c_nn!}T_n(\xi /\sqrt{2})$ where $c_0=1$ and $c_n=2$ for $n\ge 1$. Since Chebyshev polynomials are essentially cosine functions,

  $$\begin{aligned} T_j(x)T_l(x)=\frac{1}{2}T_{j+l}(x)+\frac{1}{2}T_{\vert j-l\vert }(x). \end{aligned}$$

  (A.7)

  Thus

  $$\begin{aligned} b(i,j,l)=\left\{ \begin{array}{ll} 1 &{}\quad \text {if }i=j,l=0\text { or }i=l,j=0\\ \frac{1}{2}\sqrt{\frac{c_jc_l}{c_i}}\sqrt{\frac{j!l!}{i!}} &{}\quad \text {if }j,l>0\text { and }i=j+l\text { or }i=\vert j-l\vert \\ 0 &{}\quad \text {otherwise}\\ \end{array}\right. . \end{aligned}$$

  (A.8)

  Here the expansion coefficients have a sparse pattern. For fixed j and l, there are at most two values of i such that b(i, j, l) is nonzero.

In general, we compute b(i, j, l) by matching the monomial coefficients on the both sides of (A.2) (see e.g., [38]). Suppose that

$$\begin{aligned} \varphi _n(\xi )=\sum _{m=0}^nP_{m,n}\xi ^m. \end{aligned}$$

According to (A.2), for $i>j+l$, $b(i,j,l)=0$, and $\{b(i,j,l):0\le i\le j+l\}$ satisfies the following linear system

$$\begin{aligned} \sum _{i=0}^{j+l}b(i,j,l)P_{m,i}=\sum _{r=\max \{0,i-l\}}^{\min \{i,j\}}P_{r,j}P_{i-r,l}. \end{aligned}$$

(A.9)

It is easy to solve (A.9) directly as $\{P_{m,n}\}_{m,n=0}^{j+l}$ is a upper triangular matrix. This procedure is applicable to any set of orthogonal polynomials.