[Paper Review] A multi-modes Monte Carlo finite element method for elliptic partial differential equations with random coefficients

Xiaobing Feng, Junxia Lin+1 more|arXiv (Cornell University)|Mar 29, 2016

Probabilistic and Robust Engineering Design12 references3 citations

TL;DR

This paper proposes a multi-modes Monte Carlo finite element method for solving elliptic PDEs with random coefficients that are small perturbations of deterministic coefficients. By expanding the solution in a power series of the perturbation parameter and leveraging repeated LU decompositions for identical deterministic operators across modes, the method achieves computational complexity comparable to solving only a few deterministic problems, significantly improving efficiency over standard Monte Carlo or stochastic Galerkin methods.

ABSTRACT

This paper develops and analyzes an efficient numerical method for solving elliptic partial differential equations, where the diffusion coefficients are random perturbations of deterministic diffusion coefficients. The method is based upon a multi-modes representation of the solution as a power series of the perturbation parameter, and the Monte Carlo technique for sampling the probability space. One key feature of the proposed method is that the governing equations for all the expanded mode functions share the same deterministic diffusion coefficients, thus an efficient direct solver by repeated use of the $LU$ decomposition matrices can be employed for solving the finite element discretized linear systems. It is shown that the computational complexity of the whole algorithm is comparable to that of solving a few deterministic elliptic partial differential equations using the $LU$ director solver. Error estimates are derived for the method, and numerical experiments are provided to test the efficiency of the algorithm and validate the theoretical results.

Motivation & Objective

• To address the high computational cost of solving elliptic PDEs with random coefficients using standard Monte Carlo or stochastic Galerkin methods.
• To develop a method that efficiently handles weakly random diffusion coefficients of the form $ a_0 + \varepsilon \eta $, where $ \varepsilon $ is small.
• To exploit the structural similarity of mode equations—sharing the same deterministic operator—to enable reuse of LU decompositions across modes.
• To achieve computational efficiency comparable to solving a few deterministic problems, even when sampling over a high-dimensional probability space.
• To extend the framework to general random media via Karhunen-Loève expansion, enabling broader applicability.

Proposed method

• The solution $ u^\varepsilon $ is expanded as a power series in $ \varepsilon $: $ u^\varepsilon = \sum_{n=0}^\infty \varepsilon^n u_n $, where each mode $ u_n $ satisfies an elliptic PDE with deterministic coefficients and a random source.
• The Monte Carlo method is applied to sample the probability space for the random field $ \eta $, while the finite element method discretizes each mode equation.
• The finite element stiffness matrix for each mode is identical due to the shared $ \nabla \cdot (a_0 \nabla) $ operator, enabling reuse of LU decomposition matrices across all modes.
• For each sampled realization $ \omega_j $, solving $ u_n(\omega_j) $ reduces to forward and backward substitutions using precomputed $ L $ and $ U $ matrices, drastically reducing solve time.
• The method is generalized to broader random media by rewriting the coefficient as $ \overline{a}(x) + \varepsilon \zeta(x,\omega) $ via Karhunen-Loève expansion, where $ \varepsilon = O(\ell^{r/2}) $.
• The approach is flexible and can be combined with any spatial discretization method, including finite differences, discontinuous Galerkin, or spectral methods.

Experimental results

Research questions

• RQ1 Can a multi-modes expansion combined with Monte Carlo sampling reduce the computational cost of solving random elliptic PDEs with weakly random coefficients?
• RQ2 How can the shared deterministic operator across all mode functions be exploited to accelerate solution of the linear systems arising from finite element discretization?
• RQ3 What is the computational complexity of the proposed method compared to classical Monte Carlo or stochastic Galerkin methods?
• RQ4 To what extent can the method be generalized to random coefficients that are not initially in the form $ a_0 + \varepsilon \eta $?
• RQ5 How does the accuracy of the approximation depend on the number of modes $ N $ and the perturbation size $ \varepsilon $?

Key findings

• The computational complexity of the proposed multi-modes Monte Carlo finite element method is comparable to solving only a few deterministic elliptic problems using an LU direct solver, due to reusable LU decompositions across modes.
• For $ \varepsilon = 0.2 $, the relative $ L^2 $-norm error drops from 0.0104 (with $ N=2 $ modes) to 0.0016 (with $ N=5 $), indicating high accuracy with moderate mode count.
• When $ \varepsilon = 0.4 $, the error decreases from 0.0416 ($ N=2 $) to 0.0016 ($ N=5 $), showing improved convergence with increasing mode count.
• For $ \varepsilon = 0.8 $, the error reduces from 0.1632 ($ N=2 $) to 0.0138 ($ N=5 $), confirming that more modes are needed for larger perturbations.
• The error platetus at $ \varepsilon = 0.2 $ for $ N=5 $ is due to spatial discretization error dominating the total error, not the mode expansion.
• The method remains effective when generalized to general random media via Karhunen-Loève expansion, by expressing the coefficient as $ \overline{a}(x) + \varepsilon \zeta(x,\omega) $ with $ \varepsilon = O(\ell^{r/2}) $, enabling application to Gaussian random fields with correlation length $ \ell $.