A Novel Method for Solving Multiscale Elliptic Problems with Randomly Perturbed Data

We propose a method for efficient solution of elliptic problems with multiscale features and randomly perturbed coefficients. We use the multiscale finite element method (MsFEM) as a starting point and derive an algorithm for solving a large number of multiscale problems in parallel. The method is intended to be used within a Monte Carlo framework where solutions corresponding to samples of the randomly perturbed data need to be computed. We show that the proposed method converges to the MsFEM solution in the limit for each individual sample of the data. We also show that the complexity of the method is proportional to one solve using MsFEM (where the fine scale is resolved) plus N (number of samples) solves of linear systems on the coarse scale, as opposed to solving N problems using MsFEM. A set of numerical examples is presented to illustrate the theoretical findings.


Introduction.
Multiscale problems are some of the greatest challenges in computational mathematics today.In all branches of the engineering sciences, we encounter problems with features on several different scales.Flow in a porous medium, such as the earth's subsurface, is a typical example.Many of these problems can be modeled using partial differential equations with multiscale features in the coefficients.In practice, these coefficients are often the result of experimental and/or field measurements, which causes some level of uncertainty.This uncertainty can be modeled as random perturbations.
We consider the problem of finding u satisfying (1.1) where Ω ⊂ R dim is a polygonal domain with boundary ∂Ω, where dim = 1, 2, 3, A(x) is a scalar stochastic function governed by some probability structure, and f (x) is a given deterministic function.In particular, we assume that samples of A are positive, bounded coefficients exhibiting heterogeneity in multiple scales with measurement errors and f ∈ L 2 (Ω).
In order to get a reliable solution to (1.1) we need to address two problems: 1.The coefficient A varies over different scales in space.Resolving the finest scale on a single mesh would yield a huge number of unknowns.This calls for multiscale methods.The idea is to divide the fine scale field into many local subproblems and solve these in order to modify a global coarse scale equation.This leads to a coarse scale equation in which the fine scales are taken into account.2. The coefficient A has measurement errors which can be modeled as a random perturbation around a mean.This yields a distribution of A as input data and a distribution of u as output data.Obtaining access to stochastic quantities of interest of the output is computationally challenging.During the last decade several multiscale methods have been developed which deal with the first difficulty, see, e.g., [3,11,19,21,24].In these works two scales are mainly considered; a coarse scale, on which the multiscale features are not resolved but where computations can be performed globally, and a fine scale, where the multiscale features are resolved but where global computations are not possible.The fine scale equations are decoupled and solved locally.The result is used to modify the coarse scale equations so that they include information from the fine scale.
The second issue has also been studied extensively.A common approach is the perturbation method.Here the random functions and operators in the differential equation are expanded in a Taylor series about their respective mean values (see [23]).Another approach is to expand the random function in a Karhunen-Loève expansion (see, e.g., [7,18]).A third approach often used in applications due to its simplicity is Monte Carlo simulation, i.e., to solve (1.1) for a large number of samples of the coefficient A in order to get the distribution of the solution u.
Recently, works that consider both issues simultaneously have been presented (see, e.g., [1,8,25]).In [1] a small number of realizations A are selected, and the multiscale finite element method (MsFEM) is used to compute corresponding multiscale basis functions.Then given an arbitrary realization A, the solution u is sought within the span of the precomputed multiscale basis functions.In this paper we attack a related issue.We consider a Monte Carlo type situation where a number of samples of A are given and we wish to compute corresponding solutions u by solving (1.1) using MsFEM.We present an efficient way of doing this for a vast number of samples by using the fact that the local matrices we need to invert on the fine scale are very similar in structure among the samples.This idea was first presented in a series of papers [15,16,17] where nonoverlapping domain decomposition was used to localize the computations.A great advantage of using MsFEM instead of domain decomposition is that no iteration is needed and no communication between the subdomains is needed.
The method depends on an assumption that the random variation of the elliptic coefficient is described as a piecewise constant (or polynomial) function on the coarse mesh.The multiscale basis functions associated with a coarse grid are represented as a truncated Neumann series expansion.We give a theoretical convergence analysis that shows that this series converges geometrically, with respect to the number of terms in the Neumann series, to the exact multiscale basis functions as we add more terms to the series.A similar convergence property holds for the global solution computed with this representation approximating the global MsFEM solution.The complexity of the method, when solving for N samples of the data A, is proportional to solving once using standard MsFEM plus solving N linear systems of equations on the coarse scale.This should be compared to solving N problems using standard MsFEM.
The rest of the paper is organized as follows.In section 2 we present modeling assumptions and the computational method.We derive a rigorous convergence analysis of the method in section 3. Section 4 gives a set of numerical examples to confirm the theoretical results.

Modeling assumption and computational method.
In this section we focus on modeling assumptions and present the method we use for computing samples of the solution to (1.1).We start by establishing some notation.
We let T h be a set of finite elements τ discretizing Ω, i.e., Ω = τ ∈T h τ , h τ being the diameter of τ and h = max τ ∈T h h τ .Each element τ has d vertices, i.e., we assume that τ are simplices.Finally, for a point in R 2 , we use the notation (x, y), which should be clearly distinguishable from the L 2 scalar product by examining the context.

Weak formulation and modeling assumption.
We formulate the variational formulation of (1.1): We assume that the stochastic elliptic coefficient can be decomposed into deterministic and random components, A = a + A. Further we assume that there exist constants 0 < a 0,τ ≤ a 1,τ < ∞ so that a 0,τ ≤ a(x) ≤ a 1,τ for all x ∈ τ for all τ ∈ T h and that −δ 0 a 0,τ ≤ A τ ≤ δ 1 a 0,τ < ∞ for all τ ∈ T h , where 0 ≤ δ 0 < 1 and 1 ≤ δ 1 < ∞.A simple application of the Lax-Milgram theorem guarantees the existence and uniqueness of the corresponding solution for each sample of A, and furthermore, for some constant c > 0 depending on a 0,τ , δ 0 , and Ω.We assume that the random perturbation A has a particular structure.Let where χ τ (x) denotes the characteristic function on τ .With this structure, the variational formulation for (1.1) is to seek u ∈ H 1 0 (Ω) satisfying where for an arbitrary domain D ⊆ Ω.

VICTOR GINTING, AXEL M ÅLQVIST, AND MICHAEL PRESHO
Within the framework of Monte Carlo simulation for quantifying uncertainty effects of this random coefficient, the basic idea relies on treating (2.4) as a black box producing u given finite samples of A(x).The realization data can then be used to calculate the probability distribution and/or some statistics of a functional (quantity of interest) of u.In practice, one employs the finite element method to find the approximation of u belonging to certain finite dimensional subspace of H 1 0 (Ω).The finite element method reads: find where V h , e.g., could be the space of continuous piecewise linear functions on T h vanishing on ∂Ω.
A heuristic consideration of accuracy suggests that the elliptic coefficient A dictates the mesh configuration T h .In this case, the deterministic coefficient a can exhibit heterogeneity in several different length scales of several orders of magnitude.Consequently, for an accurate solution, T h is constructed so that the heterogeneity of a is captured.This leads to a large linear system of equations that needs to be solved for each sample of A. Obviously this fact requires a tremendous amount of computational work when a large number of samples are used in the Monte Carlo simulation.In the next section, we discuss one alternative to overcome this, namely, the use of the MsFEM.
From now on, we only consider samples of the stochastic function A. In particular, the method we present is designed so that samples of the corresponding solutions u are computed in parallel with no interaction.This means that both the method and the analysis can be performed on a single sample of A without loss of generality.For this reason we will, from now on, refer to A as a particular realization of the data.In the same way, u will be the corresponding realization of the solution.The reason for this slight abuse of notation is that the paper will be easier to follow.

The MsFEM.
The MsFEM was introduced in [19] and further analyzed in [20].In its basic form, MsFEM is a finite element method that is based on solving the variational equation (2.6) in a finite dimensional subspace of H 1 0 (Ω).The distinction is in its aim to compute the solution on a coarse scale (posed on coarser mesh T h ) without directly resolving the fine scale heterogeneity globally.The key ingredient is the construction of an appropriate multiscale finite dimensional space in which the solution is sought.In particular, the fine scale heterogeneity in a should be imbedded in this finite dimensional space.This information is incorporated into coarse scale formulation through the global/coarse scale stiffness matrix.We point out that the idea of using basis functions satisfying certain differential equations has been used before (see, e.g., [4,5,6] and references therein).
We set the so-called multiscale basis functions φ i,τ for vertices i = 1, . . ., d of finite element τ to satisfy (2.7) where g i,τ (x) is a piecewise linear/bilinear function on ∂τ with g i,τ (x j ) = δ ij , x j being the vertex coordinate of τ .Then the multiscale finite dimensional space is defined as The MsFEM solution is u ms,h ∈ V ms,h satisfying (2.9) Application of this procedure is summarized in Algorithm 1.
Algorithm 1 MsFEM.Construct the multiscale basis functions {φ i,τ , i = 1, . . ., d} for every At this stage we define a projection operator Q ms,h : Using this projection we have the identity u ms,h = Q ms,h u.

A Neumann series approximation of MsFEM basis functions.
Ms-FEM appropriately addresses the first problem described in the introduction.However, within the framework of Monte Carlo simulation, we want to solve u associated with a given realization of A. Thus, for a large number of realizations which is common in Monte Carlo simulation, building up the statistics for u translates into many executions of Algorithm 1.In particular, we need to compute the multiscale basis functions for a given realization of A. This can easily become very expensive due to the slow convergence of the Monte Carlo algorithm.
In a series of papers [15,16,17], a technique using truncated Neumann series was introduced for Lions' nonoverlapping domain decomposition algorithm.In this section we propose that a similar idea is directly applicable for MsFEM and has the potential for an efficient and robust numerical procedure when applied to problems described in section 2. In this setting, the multiscale basis functions in MsFEM will be approximated by a truncated Neumann series.Using this representation, computation of the multiscale basis functions can be reduced to a one-time preprocessing effort.We now describe the procedure in detail.
We let Π : Notice that by letting , where I is the identity, and by using (2.7), we may write a variational formulation: seek

VICTOR GINTING, AXEL M ÅLQVIST, AND MICHAEL PRESHO
At this stage, we are in position to construct the Neumann series multiscale basis functions.To offer some motivation, we consider a linear algebra example associated with the variational formulation (2.11) which can be written as if S has an inverse.Given a constant λ τ , we multiply the last equation by 2/(2+A τ λ τ ) and rewrite it appropriately as Under appropriate assumptions the solution can be written as a Neumann series, Our intention is to apply the linear algebra exercise in formulating the Neumann series multiscale basis functions within the variational formulation (2.11).We propose Algorithm 2.

Algorithm 2 Neumann Series Multiscale Basis Functions
Fix i and τ , and set a positive constant λ τ .Compute φ Note that when λ τ = 0, this procedure is similar to the one introduced in [15] for nonoverlapping domain decomposition.We will show later that a judicious choice of λ τ will guarantee convergence of the series for any A τ fulfilling −δ 0,τ a 0,τ ≤ A τ ≤ δ 1,τ a 0,τ , where 0 ≤ δ 0 < 1 and δ 1 ≥ 1.We also note that similar expansion has been introduced in [14].Once this is in place, the corresponding multiscale finite dimensional space is defined as The MsFEM with truncated Neumann series solution is given by the following: find u ms,h ∈ V ms,h such that We Again we define a projection operator Q ms,h : We have the identity u ms,h = Q ms,h u.

Complexity.
Here we give a brief description of the amount of work required for the MsFEM Neumann series approximation.We think of the complexity in the framework of the Monte Carlo procedure, where N samples need to be computed.
We let |T h | be the number of coarse elements τ ∈ T h , n τ be the degrees of freedom on the subgrid of element τ (for simplicity we assume the same number for all coarse elements), and n T be the degrees of freedom on the coarse mesh T h .Furthermore, we let L(n) be the amount of work needed to compute the solution to one linear system with n unknowns.
We start by computing the work needed for matrix assembly from one coarse element τ .The coarse entries that need to be computed are a τ ( φi,τ , φj,τ ), b τ ( φi,τ , φj,τ ), and (f, φj,τ ) for 1 ≤ i, j ≤ d.In order to minimize the work it is crucial not to compute the basis functions { φi,τ } d i=1 first and then compute the products needed for the coarse scale computation.Instead we expand φi,τ in terms of φ(p) i,τ and (I − Π)g i,τ and compute terms of the form a ), and so on for 1 ≤ i, j ≤ d and 1 ≤ p, q ≤ P.This way we avoid multiplying the random numbers A τ with vectors of size n τ .Instead we multiply A τ with numbers.
The amount of work needed to compute φ(p) i,τ and Πg i,τ for all 1 ≤ i ≤ d and 1 ≤ p ≤ P is of the order d • P • L(n τ ).The amount of work needed to compute the products of type a τ ( φ(p) i,τ , φ(q) j,τ ) for all 1 ≤ i, j ≤ d and 1 ≤ p, q ≤ P is of the order This means that the amount of work needed before we involve the random numbers is proportional to Given the numbers a τ ( φ(p) i,τ , φ(q) j,τ ) and the corresponding numbers for b τ (•, •) and (f, •) for all 1 ≤ i, j ≤ d and 1 ≤ p, q ≤ P, we can assemble the matrices and vectors for all N samples.This work will be proportional to N is the amount of work needed to compute (−2A τ /(2 + λ τ A τ )) 2P (worst case) given a single realization A τ , i.e., it is proportional to P. If we do this for all coarse elements τ the work scales like 3 .Finally we need N global solves on the coarse grid at the cost N • L(n T ).In total we get (2.16) The with some moderate-sized constant.Equation (2.17) should be compared to solving for all samples using MsFEM without truncation which would be Note that N typically scales like one over the subgrid mesh size squared in order to balance the discretization error (proportional to the subgrid mesh size for energy norm estimates) and the statistical error (proportional to N −1/2 according to the CLT) in a computation.(See [15,16] for a more extensive discussion.)This means that N typically needs to be very large if we want to get convergence for some statistical output quantity of interest.

Piecewise polynomial perturbation.
The method can easily be extended to piecewise polynomial perturbations.This makes it possible to consider continuous perturbations which can be more natural if a is continuous in itself.We give a brief description of how the algorithm is modified when the perturbation is not piecewise constant.
On each element τ we assume that A τ = m j=1 A τ j ϕ j , where {ϕ i } m i=1 is a partition of unity that forms a basis for polynomials of a certain degree on τ .The standard MsFEM basis functions would then fulfill where b j τ (v, w) = (ϕ j ∇v, ∇w) τ .If, for simplicity, we let λ τ = 0 and study the equivalent linear algebra example, we get This series can again be truncated (in p).The computational cost increases since we need to solve P • m P + 1 systems on each element for each basis function initially, instead of P + 1.This is because the number of terms in the highest power of the sum is m P , and for each term P linear systems need to be solved.The addition 1 is in order to compute the projection Πg i,τ .The entries in the lower powers p are given during the computation of the highest power P. It is clear that m and P need to be fairly small numbers; otherwise, this procedure will be very costly.Given this formulation, an algorithm which is very similar to Algorithm 2 can be constructed.

2.7.
Oversampling.One drawback of MsFEM is the error resulting from imposing an artificial boundary condition g i,τ when computing the multiscale basis function (2.7).This creates a discrepancy which in general can be quantified by the ratio of the physical scale to the coarse mesh size.An analysis of the linear MsFEM, which shows that the convergence depends on this ratio, has been done in [20].Oversampling is a technique to overcome this drawback, in which the local problems (2.7) are solved in domains whose size is larger than τ .This procedure has been proposed and analyzed in [13] for linear MsFEM.
In the context of Neumann series approximation, oversampling yields a situation in which A τ is a piecewise constant function in the domain where the multiscale basis functions are computed.Ideally, this is taken into account by introducing more parameters in the local approximation of A τ , one for each neighboring element, in the same fashion as in the section on nonpiecewise constant perturbation above.This technique has also been studied in [17].A simpler approach is to approximate A τ by its value on the original coarse element τ .We intend to study these approaches in the future.

Error analysis.
The aim of this error analysis is primarily to show the convergence of the Neumann series and to estimate the error introduced by truncating the Neumann series in the global MsFEM solution.Thus, the expectation from this analysis is to devise a way for controlling the error produced by the Neumann series approximation in terms of the method parameters.We start by proving that the Neumann series approximation of the multiscale basis functions converges on each coarse element τ .We first prove the following three technical lemmas.
Proof.The proof is a straightforward calculus exercise which we present in the Appendix.
We now turn to the global solution.We show that u ms,h converges to u ms,h as P → ∞, and we also compute a rate for this convergence.We start with a lemma.
VICTOR GINTING, AXEL M ÅLQVIST, AND MICHAEL PRESHO Lemma 3.5.Let (A τ , λ τ ) ∈ A λ,τ , and assume that there exist constants 0 ≤ δ 0 < 1 and δ 1 ≥ 1 such that −δ 0 a 0,τ ≤ A τ ≤ δ 1 a 0,τ .Then for some coefficients α i,τ 's and β i,τ 's.Using these coefficients, we define We note that by using (2.12) and Theorem 3.4 and by denoting γ With all these identities and setting We first bound (Q ms,h u− w) and ( Q ms,h u−w).Using (3.6) and applying the Cauchy-Schwarz inequality and Lemma 3.2 gives where we have used a τ (g i,τ − Πg i,τ , z) = 0. We apply these to get the estimate where we have also used |||Q ms,h u||| ≤ |||u||| and the definition of in the statement of the lemma.Next, we take the norm of (3.7), apply Lemma 3.2, and add and subtract Since |γ| P c P λ,τ < 1, we can collect the second term on the right side to the left side of this last inequality.Moreover, using a similar argument for deriving the estimate of z τ , we deduce that Q ms,h u − which in turn yields Finally, since both Q ms,h and Q ms,h are projections, we have The conclusion of the lemma follows immediately by applying (3.8) and (3.10) to this last inequality.
We are now ready to present the main result of the analysis.Theorem 3.6.Let u ∈ H 1 0 (Ω) be the solution of (1.1), u ms,h ∈ V ms,h be the solution of (2.9), and u ms,h ∈ V ms,h be the solution of (2.14).Then Copyright © by SIAM.Unauthorized reproduction of this article is prohibited.

VICTOR GINTING, AXEL M ÅLQVIST, AND MICHAEL PRESHO
Using the definition of the projections Q ms,h and Q ms,h we write Rearranging these identities and applying the Cauchy-Schwarz inequality, we get where we use Lemma 3.5 in the last step.This concludes the proof.Corollary 3.7.The solution to the modified system u ms,h converges to u ms,h as Proof.This follows immediately from Theorem 3.6 since < 1.
We emphasize here that Theorem 3.6 shows the robustness of the Neumann series multiscale basis functions within the framework of MsFEM.This has been accomplished by comparing the usual MsFEM solution with the MsFEM solution using this truncated series.The comparison results in bounding the difference between the two solutions in terms of known constants which can also be made arbitrarily small by choosing sufficiently large P, i.e., more terms in the series.
Having said this, the fact remains that the MsFEM solution has its own error as compared to the true solution of the original problem (1.1).Several attempts have been made to derive errors in the MsFEM solution.(See [12,19,20] for details.)It suffices to say that most of the analyses apply homogenization techniques which restrict the elliptic coefficient to periodic functions.Nevertheless, MsFEM has been tested extensively during the last decade and has been shown to perform well for many applications.(See, for example, [9,10,11,22].)Below we present a theorem where we bound the total error in terms of the error in the standard MsFEM solution and in the error committed by using the truncated Neumann series.
Proof.We use the triangle inequality to get We refine the analysis from Theorem 3.6 as follows.Consider I 1 and I 2 in the proof of We combine (3.13) and (3.15) and use Lemma 3.5 to get the desired result.Remark 3.1.For simplicity we have considered P to be a fixed number.It can easily be chosen differently on different coarse elements τ based on the local size of a 1,τ /a 0,τ .This will not affect the results in the analysis, but it will make the presentation more cumbersome.
Remark 3.2.When comparing Theorems 3.6 and 3.8, it appears that the order in terms of differs.The trick from (3.14) could be used in the proof of Theorem 3.6 as well, leading to terms of the form |||u − Q ms u||| instead of |||u|||.But since we have no assumption on the quality of the standard MsFEM solution Q ms u and no possibility to find a similar term in that left-hand side to subtract against we have decided to leave it out.This means that if Q ms u is close to u, the bound in Theorem 3.6 might be too pessimistic.
Remark 3.3.The convergence of the method suffers clearly when gets close to 1, i.e., when the perturbations are large compared to the underlying deterministic function a(x).This means that the method works best if the randomness is a perturbation around some mean.In practice this corresponds to a situation where we have a fairly good understanding of the main features of the diffusion coefficient a(x) but the data are associated with measurement errors.If the measurements are very sparse and the coefficients in the computation are basically given as realizations of some stochastic model, this will not be achieved and the convergence will be slow.convergence of the multiscale basis functions.For this purpose, we use an element τ = [0, 1] × [0, 1] in which we compute a nodal basis function associated with the lower left vertex.In this figure, the left profile is a periodic coefficient such that 0.25 ≤ a(x) = (1 − 0.99 sin(10πx 1 )) −1 (1 − 0.99 sin(10πx 2 )) −1 ≤ 10 4 with an average of 62.4, and the right profile is a coefficient exhibiting channelized features with 0.24 ≤ a(x) ≤ 15459.1 and an average of 108.865.The multiscale basis functions are computed by discretizing τ into 100 × 100 squares for the periodic coefficient and into 120 × 120 squares for the channelized coefficient.A comparison is made between the usual multiscale basis functions and the truncated Neumann series multiscale basis functions.A standard linear finite element is employed to compute these functions.
As described in the proof of Theorem 3.4, the error associated with the Neumann series approximation of the multiscale basis functions is bounded by A τ [I(A τ , λ τ )] P , where One way to attempt to numerically observe the dependency of λ τ is to fix A τ and vary λ τ accordingly.
We note that the deterministic periodic coefficient has 1/a 0,τ + 1/a 1,τ ≈ 3.9, so by consulting Lemma 3.3, we fix A τ = 2 and vary λ τ ≥ 3.9.The left plot of Figure 4.2 shows the values of I(2, λ τ ).The middle and right plots show the error as a function of I(2, λ τ ) for P = 4 and P = 9, respectively, in log-scale.A linear regression done for these data reveals that the slopes are 4.8 and 10.7.This confirms the theoretical convergence in Theorem 3.4.Figure 4.3 shows a similar assessment for a deterministic channelized coefficient.Here 1/a 0,τ +1/a 1,τ ≈ 4.1, so we fix A τ = 1 and vary λ τ ≥ 4.1.deterministic coefficients.All plots are in logarithmic scale, and the slopes in each of them are obtained by linear regression of the data involved.Again, the slopes in these figures confirm the theoretical finding in Theorem 3.4.

Convergence of MsFEM solutions using Neumann series basis functions.
In this subsection we present a numerical experiment to show convergence of MsFEM using the Neumann series basis functions.Here we compare u h (MsFEM solution using the usual multiscale basis functions) and u h (MsFEM solution using the Neumann series basis functions).We solve the elliptic problem   The averages of a(x) in the coarse elements range from 0.44869 to 181.359 for the periodic coefficient and from 4.39 to 1114.02 for the channelized coefficient.Perturbation in each element is 30% of its average for the periodic coefficient and 50% of its average for the channelized coefficient.The plots indicated by "Reference" are the MsFEM solutions using the usual multiscale basis functions.We can see that as we add more terms in the series, the approximate solutions do converge to the reference solutions.

Conclusion.
In this paper, we have proposed an efficient procedure for solving multiscale elliptic problems with randomly perturbed coefficients.The theoretical perspective that we offer confirms that the Neumann series approach for computing the multiscale basis functions in MsFEM converges geometrically with respect to the number of terms in the series.Similarly, the MsFEM solution computed using these Neumann series converges to the MsFEM solution computed using the usual multiscale basis functions.The complexity of the method when solving N samples of the data is proportional to using standard MsFEM on a single sample plus solving N linear systems of equations on the coarse scale.
We wish to investigate a number of open problems in the future.As mentioned in section 2, an interesting extension to this method is to use the oversampling in computing the multiscale basis functions.Moreover, we would like to apply the method to several applications, e.g., in oil reservoir simulation.One idea is to quantify uncertainty in the carbon production due to the randomness of the measured data of the reservoir's permeability.This would be done within the Monte Carlo framework, and we aim at including statistical errors for a given output quantity in the analysis.

4 .
Numerical examples.We present several numerical experiments which illustrate the theoretical investigations described in the previous sections.Our examples are mainly used to study the convergence of the Neumann series in MsFEM.4.1.Convergence of the Neumann series multiscale basis functions.

Figure 4 .
1 shows two deterministic elliptic coefficients a(x) that we use to test the

Fig. 4 . 6 .Fig. 4 . 7 .
Fig. 4.6.Behavior of the error of the MsFEM solution using the Neumann series multiscale basis function with respect to for periodic (left) and channelized (right) deterministic elliptic coefficients.
Figure 4.6 shows the convergence for the two deterministic coefficients.We use P = 4 for this example.The asymptotic slopes for both results are 5.3 and 5.5, respectively.

Figure 4 .
7 shows a comparison of the MsFEM solutions for several numbers of terms in the Neumann series, using periodic and channelized deterministic coefficients.
may now modify Algorithm 1 to include the Neumann series approximation (see Algorithm 3).Algorithm 3 MsFEM with Truncated Neumann Series For every τ ∈ T h and i = 1, . . ., d call Algorithm 2 to get φ Construct the multiscale basis functions { φ i,τ , i = 1, . . ., d} for every τ ∈ T h with constants we have neglected are independent of the other quantities in the calculation d, P, n τ , n T , |T h |, and N .For small values of P (which is what we typically are interested in) and small values of d (typically 3 or 4), we essentially get that the work is proportional to(2.17) Theorem 3.6 and estimate using the Cauchy-Scwharz inequality, and |ab| ≤ β 2 |a| 2 + h u||| 2 from both sides, divide with (1− β 2 ), take square roots, and use |x| 2 + |y| 2 ≤ |x| + |y| to get