# Using subspace substitution to obtain rapidly convergent series expansions
for a class of resolvents

Graeme W. Milton Department of Mathematics, University of Utah, e-mail:
<EMAIL_ADDRESS>
###### Abstract

Following advances in the abstract theory of composites, we develop rapidly
converging series expansions about $z=\infty$ for the resolvent ${\bf
R}(z)=[z{\bf I}-{\bf P}^{\dagger}{\bf Q}{\bf P}]^{-1}$ where ${\bf Q}$ is an
orthogonal projection and ${\bf P}$ is such that ${\bf P}{\bf P}^{\dagger}$ is
an orthogonal projection. It is assumed that the spectrum of ${\bf
P}^{\dagger}{\bf Q}{\bf P}$ lies within the interval $[z^{-},z^{+}]$ for some
known $z^{+}\leq 1$ and $z^{-}\geq 0$ and that the actions of the projections
${\bf Q}$ and ${\bf P}{\bf P}^{\dagger}$ are easy to compute. The series
converges in the entire $z$-plane excluding the cut $[z^{-},z^{+}]$. It is
obtained using subspace substitution, where the desired resolvent is tied to a
resolvent in a larger space and ${\bf Q}$ gets replaced by a projection
$\underline{{\bf Q}}$ that is no longer orthogonal. When $z$ is real the rate
of convergence of the new method matches that of the conjugate gradient
method.

Keywords: Resolvent, Series Expansions, Matrix Inverses, Subspace
Substitution, Iterative Methods.

MSC codes: 65B10, 65F10, 65Z05, 47A10, 35B27

##  1 Introduction

Here we seek rapidly converging series expansions about $z=\infty$ of
resolvents of the form

${\bf R}(z)=[z{\bf I}-{\bf P}^{\dagger}{\bf Q}{\bf P}]^{-1},$ (1.1)

where

  * •

${\bf Q}:{\mbox{\boldmath${\cal H}$}}\to{\mbox{\boldmath${\cal Q}$}}$ is an
orthogonal projection from a Hilbert or vector space ${\cal H}$ onto
${\mbox{\boldmath${\cal Q}$}}\subset{\mbox{\boldmath${\cal H}$}}$,

  * •

For some subspace ${\mbox{\boldmath${\cal E}$}}\subset{\mbox{\boldmath${\cal
H}$}}$, ${\bf P}$ (with adjoint ${\bf P}^{\dagger}$) is a mapping from ${\cal
E}$ to ${\cal H}$ (and vice-versa),

  * •

${\bf P}$ is such that $\mbox{\boldmath${\Gamma}$}={\bf P}{\bf
P}^{\dagger}:{\mbox{\boldmath${\cal H}$}}\to{\mbox{\boldmath${\cal H}$}}$ is
an orthogonal projection onto ${\cal E}$,

  * •

${\bf P}$ is such that ${\bf P}^{\dagger}{\bf P}={\bf
I}:{\mbox{\boldmath${\cal E}$}}\to{\mbox{\boldmath${\cal E}$}}$ is the
identity on ${\cal E}$.

To begin, we keep things general and consider the resolvent

${\bf R}(z)=[z{\bf I}-{\bf P}^{\dagger}{\bf A}{\bf P}]^{-1},$ (1.2)

where ${\bf A}:{\mbox{\boldmath${\cal H}$}}\to{\mbox{\boldmath${\cal H}$}}$ is
Hermitian satisfying the bounds $a^{-}{\bf I}\leq{\bf A}\leq a^{+}{\bf I}$ for
some known real $a^{+}$ and $a^{-}$ that could be negative. Such resolvents
appear in the abstract theory of composites and their computation using series
expansions has been the subject of much attention: see [29] and references
therein. It is not the purpose of this paper to review the field. Rather the
paper aims to bring to the wider numerical analysis community an approach for
developing a rapidly converging series expansion of the resolvent in (1.1)
that has been discovered in the theory of composites. It uses the idea of
subspace substitution with non-orthogonal subspaces.

The results are based on those in Chapter 8 of [23] and [21] which in turn are
based on many works that we will cite as we proceed.

By rescaling and shifting ${\bf A}$ and $z$, and rescaling ${\bf R}$, we can
assume without loss of generality that

${\bf I}\geq{\bf A}\geq 0.$ (1.3)

Specifically, with the replacements

${\bf A}\rightarrow(a^{+}-a^{-}){\bf A}+a^{-}{\bf I},\quad
z\rightarrow(a^{+}-a^{-})z+a^{-},\quad{\bf R}\rightarrow{\bf
R}/(a^{+}-a^{-}),$ (1.4)

(1.1) is still satisfied with ${\bf A}$ now satisfying (1.3). We make the
further assumption that the resolvent ${\bf R}_{\bf A}(a)=(a{\bf I}-{\bf
A})^{-1}$ can easily be computed when $a\notin[0,1]$. In addition, we suppose
that constants $z^{+}\leq 1$ and $z^{-}\geq 0$ are known such that

$z^{+}{\bf I}\geq{\bf P}^{\dagger}{\bf A}{\bf P}\geq z^{-}{\bf I}.$ (1.5)

Multiplying (1.5) on the left by ${\bf P}$ and on the right by ${\bf
P}^{\dagger}$ gives

$z^{+}\mbox{\boldmath${\Gamma}$}\geq\mbox{\boldmath${\Gamma}$}{\bf
A}\mbox{\boldmath${\Gamma}$}\geq z^{-}\mbox{\boldmath${\Gamma}$}.$ (1.6)

In the theory of composites the case where ${\bf A}={\bf Q}$ is a projection
corresponds to having a two-phase medium. Then there is a simple formula for
the resolvent ${\bf R}_{\bf A}(a)$:

${\bf R}_{\bf A}(a)=(a{\bf I}-{\bf Q})^{-1}={\bf Q}/(a-1)+({\bf I}-{\bf
Q})/a.$ (1.7)

While any Hermitian operator ${\bf B}$ with ${\bf I}\geq{\bf B}\geq 0$ has the
factorization ${\bf B}={\bf P}^{\dagger}{\bf Q}{\bf P}$ with

${\bf Q}=\begin{pmatrix}{\bf B}&{\bf B}^{1/2}({\bf I}-{\bf B})^{1/2}\cr{\bf
B}^{1/2}({\bf I}-{\bf B})^{1/2}&{\bf I}-{\bf B}\end{pmatrix},\quad{\bf
P}=\begin{pmatrix}{\bf I}\cr 0\end{pmatrix},$ (1.8)

the difficulty is in calculating the square roots of ${\bf B}$ and ${\bf
I}-{\bf B}$. Therefore we will assume that it is ${\bf Q}$ and ${\bf P}$ (or
${\Gamma}$) which are given, not ${\bf B}$.

The resolvent naturally appears in the solution of the equation

$[z{\bf I}-{\bf P}^{\dagger}{\bf A}{\bf P}]{\bf y}={\bf b}.$ (1.9)

Multiplying on the left by ${\bf P}/z$ and noting that ${\bf P}{\bf
y}=\mbox{\boldmath${\Gamma}$}{\bf P}{\bf y}$, we get the equivalent equations

$[\mbox{\boldmath${\Gamma}$}({\bf I}-{\bf A}/z)]{\bf E}={\bf h}\text{ with
}{\bf E}={\bf P}{\bf y}=\mbox{\boldmath${\Gamma}$}{\bf E},\quad{\bf h}={\bf
P}{\bf b}/z=\mbox{\boldmath${\Gamma}$}{\bf h}.$ (1.10)

Of course, in practice one might only be interested in solving this for one
value of $z$, which then can be absorbed into ${\bf A}$, i.e. making the
replacement ${\bf A}\rightarrow{\bf A}+z{\bf I}$. However, keeping $z$ allows
us to assess the rates of convergence of different expansions. So the
equations become simply

$\mbox{\boldmath${\Gamma}$}{\bf L}{\bf E}={\bf
h},\quad\mbox{\boldmath${\Gamma}$}{\bf E}={\bf
E},\quad\mbox{\boldmath${\Gamma}$}{\bf h}={\bf h},$ (1.11)

with

${\bf L}={\bf I}+(\sigma-1){\bf A},\quad\text{where }\sigma=(z-1)/z.$ (1.12)

In a two phase conducting composite ${\bf A}$ is the projection ${\bf
Q}=\chi{\bf I}$, where $\chi({\bf x})$ is 1 in phase one and 0 in phase two.
Then ${\bf L}=\sigma\chi{\bf I}+(1-\chi){\bf I}$ is the conductivity which
takes the value $\sigma$ in phase one and 1 in phase two. The operator
${\Gamma}$ is the projection onto gradients of potentials. We will formulate
this example more precisely in the next section.

It is helpful to rewrite (1.11) as

$[c{\bf I}+\mbox{\boldmath${\Gamma}$}({\bf L}-c{\bf I})]{\bf E}={\bf h},$
(1.13)

where $c$ is a possibly complex constant. As we will see, the choice of $c$ is
typically based on bounds on the norm of certain operators. This might not be
the best choice, particularly if the bounds are not sharp. An adaptive method
of choosing $c$ has been discovered and implemented [27].

In the theory of composites $c$ is associated with the so called reference
material ${\bf L}_{0}=c{\bf I}$. Most of the analysis here extends to the case
where ${\bf L}_{0}$ is not a multiple of the identity but we take ${\bf
L}_{0}=c{\bf I}$ to simplify the analysis.

Applying ${\bf P}$ to both sides of (1.13) and recalling that ${\bf
P}^{\dagger}{\bf P}={\bf I}$ we see that

$[z{\bf I}-{\bf P}^{\dagger}{\bf A}{\bf P}]^{-1}={\bf P}[c{\bf
I}+\mbox{\boldmath${\Gamma}$}({\bf L}-c{\bf I})]^{-1}{\bf P}^{\dagger}/z.$
(1.14)

Thus, we are left with evaluating the inverse of $[c{\bf
I}+\mbox{\boldmath${\Gamma}$}({\bf L}-c{\bf I})]$ provided $c$ is such that
this inverse exists.

To define what we mean by a rate of convergence consider a series in powers of
$\alpha$:

$\mathbb{S}=\sum_{n=0}^{\infty}\alpha^{n}{\bf C}_{n}.$ (1.15)

Let us introduce the norm $\|{\bf O}\|$ of an operator ${\bf O}$ (that need
not be self-adjoint) defined by

$\|{\bf O}\|=\max_{{\bf a},\,|{\bf a}|=1}|{\bf O}{\bf a}|,$ (1.16)

where $|{\bf b}|$ denotes the norm of a vector ${\bf b}$ in our Hilbert space.
Then an estimate of the error made in truncating the series (1.15) at $n=m$ is

$S_{m}=\|\mathbb{S}-\sum_{n=0}^{m-1}\alpha^{n}{\bf
C}_{n}\|=\|\sum_{n=m}^{\infty}\alpha^{n}{\bf
C}_{n}\|\leq\sum_{n=m}^{\infty}|\alpha^{n}|\|{\bf C}_{n}\|.$ (1.17)

It is the asymptotics of this latter series that determine the rate of
convergence $\mu$. The radius of convergence of the series is

$R=\frac{1}{\limsup_{m\to\infty}\|{\bf C}_{m}\|^{1/m}}$ (1.18)

and the series will diverge if $|\alpha|>R$ and converge if $|\alpha|<R$. We
define

$\mu=\frac{|\alpha|}{R}={|\alpha|\limsup_{m\to\infty}\|{\bf C}_{m}\|^{1/m}}$
(1.19)

as the rate of convergence. Thus, if $\mu<1$ the series converges and will be
fastest when $\mu$ is small. The series will not converge if $\mu>1$.

For the various series studied here ${\bf C}_{n}={\bf C}^{n}$ for some
operator ${\bf C}$ and the rate of convergence is simply $\mu=|\alpha|\|{\bf
C}\|$. Then, of course, there is an iterative procedure to calculate. ${\bf
e}=\mathbb{S}{\bf h}$. Clearly ${\bf e}_{m}=\sum_{n=0}^{m-1}\alpha^{n}{\bf
C}^{n}{\bf h}$ satisfies the iterative relation

${\bf e}_{1}={\bf h},\quad{\bf e}_{m+1}=\alpha{\bf C}{\bf e}_{m}+{\bf
h},\,\,m=1,2,3,\ldots,$ (1.20)

and if the iterates converge, they converge to ${\bf e}$.

##  2 A class of examples that includes the conductivity and Schrödinger
equations

Equations of the form (1.11) are prevalent throughout physics. In many cases
one is interested in solving, in the entire domain of $\mathbb{R}^{3}$ or in a
unit cell of periodicity, the equation

${\bf T}(\nabla)^{\dagger}{\bf L}({\bf x}){\bf
T}(\nabla)\mbox{\boldmath${\Psi}$}({\bf x})={\bf g}({\bf x}),$ (2.21)

where ${\bf T}$ has a polynomial dependence on $\nabla$, and ${\bf g}$ is the
source term. That a host of equations in physics can be reduced to this for
was recognized by Strang [30]. While ${\bf T}(\nabla)$ is generally linear in
$\nabla$, for some equations, such as beam and plate equations, ${\bf
T}(\nabla)$ is quadratic in $\nabla$. Following [16] let ${\bf E}={\bf
T}(\nabla)\mbox{\boldmath${\Psi}$}({\bf x})$ and choose ${\bf h}$ such that
${\bf g}={\bf T}(\nabla)^{\dagger}{\bf h}$, the latter being best done in
Fourier space. Assuming ${\bf F}(\nabla)={\bf T}(\nabla)^{\dagger}{\bf
T}(\nabla)$ has an inverse, easiest to compute in Fourier space, we introduce
the projection

$\mbox{\boldmath${\Gamma}$}(\nabla)={\bf T}(\nabla)({\bf F}(\nabla))^{-1}{\bf
T}(\nabla)^{\dagger},$ (2.22)

and then (1.11) is satisfied with ${\bf h}={\bf T}(\nabla)^{\dagger}({\bf
F}(\nabla))^{-1}$ and we can recover ${\Psi}$ using

$\mbox{\boldmath${\Psi}$}=({\bf F}(\nabla))^{-1}{\bf T}(\nabla)^{\dagger}{\bf
E}.$ (2.23)

The operators $({\bf F}(\nabla))^{-1}$ and hence
$\mbox{\boldmath${\Gamma}$}(\nabla)$ are non-local in real space. Their action
is best evaluated in Fourier space where $\nabla$ is replaced by $i{\bf k}$,
where ${\bf k}$ is the wave vector. In the context of these equations bounds
$z^{+}$ and $z^{-}$ on the spectrum of ${\bf A}$ can be obtained by using, as
appropriate, quasi-convex functions [4], A-quasi-convex functions [1], or
their generalization of $Q^{*}$-convex operators [20].

The prototypical example is the conductivity equation:

${\bf j}^{\prime}({\bf x})=\mbox{\boldmath${\sigma}$}({\bf x}){\bf e}({\bf
x})-{\bf s}({\bf x}),\quad\mbox{\boldmath${\Gamma}$}{\bf e}={\bf
e},\quad\mbox{\boldmath${\Gamma}$}{\bf
j}^{\prime}=0,\quad\text{with}\quad\mbox{\boldmath${\Gamma}$}=\nabla(\nabla^{2})^{-1}\nabla\cdot,$
(2.24)

where $\mbox{\boldmath${\sigma}$}({\bf x})$ is the conductivity tensor, while
$\nabla\cdot{\bf s}$, ${\bf j}={\bf j}^{\prime}+{\bf s}$, and ${\bf e}$ are
the current source, current, and electric field, and $(\nabla^{2})^{-1}$ is
the inverse Laplacian (there is obviously considerable flexibility in the
choice of ${\bf s}({\bf x})$, the only constraints being square integrability
and that $\nabla\cdot{\bf s}$ equals the current source). An interesting twist
is that we write the equations in this form, rather than in the more
conventional form involving ${\bf j}$ directly. This is exactly what we need
to keep the left hand side of the constitutive law divergence free. As current
is conserved, $\nabla\cdot{\bf j}=\nabla\cdot{\bf s}$, implying
$\nabla\cdot{\bf j}^{\prime}=0$, which is clearly equivalent to
$\mbox{\boldmath${\Gamma}$}{\bf j}^{\prime}=0$. To show that ${\bf
e}=\nabla(\nabla^{2})^{-1}\nabla\cdot{\bf e}$, we let $V$ be the solution of
Poisson’s equation $\nabla^{2}V=-\nabla\cdot{\bf e}$ (with $V({\bf x})\to 0$
as ${\bf x}\to\infty$), i.e. $V=-(\nabla^{2})^{-1}\nabla\cdot{\bf e}$. Then
substitution gives ${\bf e}=-\nabla V=\nabla(\nabla^{2})^{-1}\nabla\cdot{\bf
e}=\mbox{\boldmath${\Gamma}$}{\bf e}$. These steps are much easier done in
Fourier space, where $\mbox{\boldmath${\Gamma}$}({\bf k})={\bf k}({\bf
k}\cdot{\bf k})^{-1}{\bf k}^{T}$. Upon identifying ${\bf L}$ with ${\sigma}$
and ${\bf E}$ with ${\bf e}$, the application of ${\Gamma}$ to the first
equation in (2.24) gives (1.11) with ${\bf h}=\mbox{\boldmath${\Gamma}$}{\bf
s}$. For a two phase conducting composite with isotropic (scalar)
conductivities we can rescale, without loss of generality, the conductivities
so phase 2 has conductivity 1. Then ${\bf L}={\bf I}+(\sigma-1){\bf Q}$ where
$\sigma$ is the conductivity of phase 1 and ${\bf Q}=\chi{\bf I}$ where $\chi$
is the characteristic function taking the value 1 in phase 1 and 0 in phase 2.

Another example is the Schrödinger equation. With $m$ denoting the mass of the
electron and $\hbar$ denoting Planck’s constant $h$ divided by $2\pi$, the
time independent Schrödinger equation for the wave function $\psi$ of a single
electron in an infinite domain with source term $g({\bf x})$, potential
$V({\bf x})$ and energy $E_{0}$ takes the form

$[V({\bf x})-E_{0}]\psi(x)-\frac{\hbar}{2m}\nabla^{2}\psi=g({\bf x}).$ (2.25)

It can be re-expressed as

$\underbrace{\begin{pmatrix}-\nabla&1\end{pmatrix}}_{{\bf
T}(\nabla)^{\dagger}}\underbrace{\begin{pmatrix}\hbar{\bf I}/2m&0\cr
0&V-E_{0}\end{pmatrix}}_{\bf L}\underbrace{\begin{pmatrix}\nabla\cr
1\end{pmatrix}}_{{\bf T}(\nabla)}\mbox{\boldmath${\psi}$}=g,$ (2.26)

and in Fourier space

$\mbox{\boldmath${\Gamma}$}=\frac{1}{1+|{\bf k}|^{2}}\begin{pmatrix}{\bf
k}\otimes{\bf k}&i{\bf k}\cr-i{\bf k}&1\end{pmatrix},\quad{\bf
E}=\begin{pmatrix}i{\bf k}\cr 1\end{pmatrix}\mbox{\boldmath${\psi}$},\quad{\bf
h}=\frac{1}{1+|{\bf k}|^{2}}\begin{pmatrix}i{\bf k}\cr 1\end{pmatrix}g.$
(2.27)

Note that ${\bf L}$ will be bounded only if $V({\bf x})$ is bounded and
coercive if only if $V({\bf x})-E_{0}\geq\epsilon$ for some $\epsilon>0$.

Now suppose $V({\bf x})$ only takes two values,

$V({\bf x})=V_{0}+\chi({\bf x})V_{1},$ (2.28)

where $\chi({\bf x})$ takes either the value $1$ or zero and $V_{0}-E_{0}$ and
$V_{0}+V_{1}-E_{0}$ are both positive ($V_{1}$ could be negative). We rewrite
(2.25) as

$\underbrace{\begin{pmatrix}-\nabla\sqrt{\hbar/2m}&\sqrt{V_{0}-E_{0}}\end{pmatrix}}_{{\bf
T}(\nabla)^{\dagger}}\underbrace{\begin{pmatrix}{\bf I}&0\cr 0&1+\chi
V_{1}/(V_{0}-E_{0})\end{pmatrix}}_{{\bf
L}}\underbrace{\begin{pmatrix}\nabla\sqrt{\hbar/2m}\\\
\sqrt{V_{0}-E_{0}}\end{pmatrix}}_{{\bf T}(\nabla)}\psi=g,$ (2.29)

where we have replaced the previous definitions of ${\bf T}(\nabla)$ and ${\bf
L}$. We now change the definitions of ${\Gamma}$, ${\bf E}$ , and ${\bf h}$ in
Fourier space to

${\Gamma}$ $\displaystyle=$ $\displaystyle\frac{1}{V_{0}-E_{0}+|{\bf
k}|^{2}\hbar/2m}\begin{pmatrix}({\bf k}\otimes{\bf k})\hbar/2m&i{\bf
k}\sqrt{(V_{0}-E_{0})\hbar/2m}\cr-i{\bf
k}\sqrt{(V_{0}-E_{0})\hbar/2m}&V_{0}-E_{0}\end{pmatrix},$ $\displaystyle{\bf
E}$ $\displaystyle=$ $\displaystyle\begin{pmatrix}\nabla\sqrt{\hbar/2m}\\\
\sqrt{V_{0}-E_{0}}\end{pmatrix}\psi,\quad{\bf h}=\frac{1}{V_{0}-E_{0}+|{\bf
k}|^{2}\hbar/2m}\begin{pmatrix}\nabla\sqrt{\hbar/2m}\\\
\sqrt{V_{0}-E_{0}}\end{pmatrix}g,$ (2.30)

where ${\Gamma}$ projects onto the subspace of fields having the same form as
${\bf E}$ as $\psi$ varies.

Then the equations still take the form (1.11) where now

${\bf L}=[{\bf I}-{\bf Q}/z],\text{ with }{\bf Q}=\begin{pmatrix}0&0\cr
0&\chi\end{pmatrix},\quad\quad z=(E_{0}-V_{0})/V_{1}.$ (2.31)

A plethora of other equations that can be expressed in the form (1.11) are
presented in [16, 17, 18, 19].

##  3 Some well known series expansions

The simplest series expansion of the resolvent (1.1) is of course

${\bf R}=\sum_{n=0}^{\infty}({\bf P}^{\dagger}{\bf A}{\bf
P})^{n}/z^{n+1}=\sum_{n=0}^{\infty}{\bf P}^{\dagger}({\bf
A}\mbox{\boldmath${\Gamma}$})^{n}{\bf P}/z^{n+1},$ (3.1)

where we have used the fact that ${\bf P}{\bf
y}=\mbox{\boldmath${\Gamma}$}{\bf P}{\bf y}$. As is well known one can
estimate the error in taking just $m$ terms in the series expansion

$\displaystyle\|{\bf R}-\sum_{n=0}^{m-1}({\bf P}^{\dagger}{\bf A}{\bf
P})^{n}/z^{n+1}\|$ $\displaystyle=$ $\displaystyle\|\sum_{n=m}^{\infty}({\bf
P}^{\dagger}{\bf A}{\bf P})^{n}/z^{n+1}\|$ (3.2) $\displaystyle\leq$
$\displaystyle\sum_{n=m}^{\infty}\|({\bf P}^{\dagger}{\bf A}{\bf
P})\|^{n}/|z|^{n+1}$ $\displaystyle\leq$
$\displaystyle\sum_{n=m}^{\infty}|z^{+}/z|^{n}/|z|=|z^{+}/z|^{m}/(|z|-z^{+}).$

Thus, the rate of convergence for the series in powers of $1/z$ is at least
$\mu_{0}=z^{+}/z$.

An improvement is the Richardson method which expands the inverse on the right
side of (1.14) giving

$\displaystyle[c{\bf I}+\mbox{\boldmath${\Gamma}$}({\bf L}-c{\bf I})]^{-1}$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}[(\mbox{\boldmath${\Gamma}$}({\bf L}/c-{\bf
I})]^{n}/c=\sum_{n=0}^{\infty}[(\mbox{\boldmath${\Gamma}$}({\bf L}/c-{\bf
I})\mbox{\boldmath${\Gamma}$}]^{n-1}({\bf L}/c-{\bf I})/c.$ (3.3)

Such series [3, 11] have played an important role in the theory of composites.
In that context ${\Gamma}$ imposes the differential constraints of the fields
and acts locally in Fourier space while ${\bf L}$ represents the material
tensor field entering the constitutive relation and it acts locally in real
space. Particularly important was the recognition by Moulinec and Suquet [26]
that fast Fourier transforms could be used to bounce back and forth between
real space and Fourier space when computing by iteration the action of the
series in (3.3) on a field. Their numerical method triggered a multitude of
further developments, see [29] and references therein, including the extension
to non-linear equations in composite materials as reviewed in [28]. Observe
that

$\|{\bf L}/c-{\bf I}\|\leq\max\\{|\sigma-c|/|c|,|1-c|/|c|\\}.$ (3.4)

The right side is minimized when $c=(1+\sigma)/2$ , which is the choice made
by Moulinec and Suquet [26]. In this case, the series expansion becomes

$[c{\bf I}+\mbox{\boldmath${\Gamma}$}({\bf L}-c{\bf
I})]^{-1}=2\sum_{n=0}^{\infty}s^{n}[\mbox{\boldmath${\Gamma}$}(2{\bf A}-{\bf
I})]^{n}/(1+\sigma)\text{ with }s=\frac{\sigma-1}{\sigma+1}.$ (3.5)

Then, since $\|2{\bf A}-{\bf A}\|\leq 1$, the rate of convergence is at least
$\mu_{1}=|s|$. If we only take $m$ terms in the series then the error can be
bounded:

$\displaystyle\|[c{\bf I}+\mbox{\boldmath${\Gamma}$}({\bf L}-c{\bf
I})]^{-1}-2\sum_{n=0}^{m-1}s^{n}[\mbox{\boldmath${\Gamma}$}(2{\bf A}-{\bf
A})]^{n}/(1+\sigma)\|$ $\displaystyle=$
$\displaystyle\|2\sum_{n=m}^{\infty}s^{n}[\mbox{\boldmath${\Gamma}$}(2{\bf
A}-{\bf I})]^{n}/(1+\sigma)\|$ (3.6) $\displaystyle\leq$ $\displaystyle
2\sum_{n=m}^{\infty}|s|^{n}\left(\|\mbox{\boldmath${\Gamma}$}\|\|2{\bf A}-{\bf
I}\|\right)^{n}/|1+\sigma|\|$ $\displaystyle\leq$ $\displaystyle
2\sum_{n=m}^{\infty}|s|^{n}/|1+\sigma|$ $\displaystyle=$ $\displaystyle
2|s|^{m}/[(1-|s|)|1+\sigma|].$

Although it will not concern us here, the rate of convergence of this series
with this value of $c$ is better than that indicated by the bound (3.6) when
$z^{-}\neq 0$. To see this, (1.6) implies

$|s|\|\mbox{\boldmath${\Gamma}$}(2{\bf A}-{\bf
I})\mbox{\boldmath${\Gamma}$}]\|\leq|s|\max\\{|2z^{-}-1|,|2z^{+}-1|\\}$ (3.7)

and thus the rate of convergence is bounded above by the right hand side of
(3.7).

If both $z^{+}$ and $z^{-}$ are known, further improvements can be obtained by
adjusting $c$ to minimize the norm of the operator
$\mbox{\boldmath${\Gamma}$}({\bf L}/c-{\bf I})\mbox{\boldmath${\Gamma}$}$
appearing on the right side of (3.3) results in a series with an even better
bound on the rate of convergence. Specifically, using (1.6) to bound this norm
gives

$\|\mbox{\boldmath${\Gamma}$}({\bf L}/c-{\bf
I})\mbox{\boldmath${\Gamma}$}\|=\|(1-c)\mbox{\boldmath${\Gamma}$}-\mbox{\boldmath${\Gamma}$}{\bf
A}\mbox{\boldmath${\Gamma}$}/z\|/|c|\leq\max\\{|(1-c)-z^{+}/z|/|c|,|(1-c)-z^{-}/z|/|c|\\}.$
(3.8)

The maximum occurs when both expressions in the maximum are equal, i.e. for

$c=1-(z^{+}+z^{-})/(2z),$ (3.9)

and in this case we get

$\|\mbox{\boldmath${\Gamma}$}({\bf L}/c-{\bf
I})\mbox{\boldmath${\Gamma}$}\|\leq\frac{z^{+}-z^{-}}{|2z-z^{-}-z^{+}|}.$
(3.10)

So the rate of convergence is at least

$\mu_{2}=(z^{+}-z^{-})/|2z-z^{+}-z^{+}|.$ (3.11)

To compare this with the convergence rate bound $\mu_{1}=|s|$ we first note
that the first series will converge only if $|s|<1$. Writing $\mu_{2}$ in
terms of $s$ gives

$\mu_{2}=\frac{(z^{+}-z^{-})|s|}{|s(1-z^{+}-z^{+})-1|}.$ (3.12)

As $s$ varies in the unit disk the minimum of $|s(1-z^{+}-z^{+})-1|$ occurs at
$s=1$ and so we get

$\mu_{2}\leq\frac{|s|(z^{+}-z^{-})}{z^{+}+z^{-}}\leq\mu_{1}.$ (3.13)

##  4 Motivating the strategy

Suppose we wanted to find an expansion that converges as rapidly as possible
of $f(z)$ about $z=\infty$ given that $f(z)$ has singularities distributed all
along the interval $[z^{-},z^{+}]$. The strategy is to map the region of the
complex plane outside this interval to the strict interior of the unit disk.
Thus we want the image of $[z^{-},z^{+}]$ to wrap entirely around the boundary
of the unit disk. First we make the fractional linear transformation from the
complex $z$-plane to the complex $\underline{\sigma}$-plane where

$\underline{\sigma}=\frac{z-z^{-}}{z-z^{+}}.$ (4.1)

This maps $[z^{-},z^{+}]$ to the entire negative real $\underline{\sigma}$
axis, and moves the expansion point from $z=\infty$ to $\underline{\sigma}=1$.
Next we introduce the mapping

$w=\sqrt{\underline{\sigma}},$ (4.2)

which takes the cut complex $\underline{\sigma}$-plane to the right half of
the $w$-plane. Now the singularities are distributed along the entire negative
imaginary axis and the expansion point remains at $w=1$. Finally, we take

$v=\frac{w-1}{w+1},$ (4.3)

which maps the right half of the $w$-plane to the unit disk in the $v$-plane
and moves the expansion point from $w=1$ to the origin in the $v$-plane.
Putting all these transformations together we see that

$v=\left(\frac{\sqrt{\frac{z-z^{-}}{z-z^{+}}}-1}{\sqrt{\frac{z-z^{-}}{z-z^{+}}}+1}\right)$
(4.4)

is the desired transformation. We will seek an expansion of the resolvent in
powers of $v$.

Putting together the inverse transformations,

$w=\frac{1-v}{v+1},\quad\underline{\sigma}=w^{2},\quad
z=\frac{z^{+}\underline{\sigma}-z^{-}}{\underline{\sigma}-1},$ (4.5)

we obtain

$z=z(v)\equiv\frac{\frac{z^{+}(1-v)^{2}}{(v+1)^{2}}-z^{-}}{\left(\frac{1-v}{v+1}\right)^{2}-1}.$
(4.6)

This may be substituted in $f(z)$ and the resulting function $f(z(v))$ can
then be expanded in powers of $v$ to obtain a rapidly converging expression
for $f$ with the expansion point at $v=0$ corresponding to the expansion point
of $f(z)$ at $z=\infty$.

##  5 Obtaining the desired series expansions

###  5.1 An accelerated series expansion that does not require knowledge of
the spectral bounds $z^{+}$ and $z^{-}$

The initial stage in obtaining a series expansion with a generally improved
rate of convergence is to obtain one where $\sqrt{\sigma}$ naturally enters.
(Later we will replace $\sigma$ by an appropriately defined
$\underline{\sigma})$). The series expansions developed in this initial stage
were introduced for the conductivity equations in composites [6] and then
extended to elasticity [13] and more general equations (Section 14.9 in [15]).
Bounds on the convergence of this accelerated scheme were derived in this
latter paper. Later it was recognized [29] that the accelerated iteration
scheme is a special case of the Douglas Rachford splitting method [5]. As
Pierre Suquet noted (private communication, July 2024), another special case
of the Douglas Rachford splitting method is the augmented Lagrangian iteration
scheme [12] previously found to have a close connection [24] with the
accelerated scheme [6]. Bounds on the convergence of this splitting method
have been derived [7].

Our starting point is again the inverse on the left side of (3.3) but now we
change the sign of $c$, which is equivalent to taking a reference tensor ${\bf
L}_{0}$ that is negative definite. Then, following [22], we have

$[-c{\bf I}+\mbox{\boldmath${\Gamma}$}({\bf L}+c{\bf I})]^{-1}=({\bf L}+c{\bf
I})^{-1}[-c({\bf L}+c{\bf I})^{-1}+\mbox{\boldmath${\Gamma}$}]^{-1}=2({\bf
L}+c{\bf I})^{-1}({\bf K}-\mbox{\boldmath${\Upsilon}$})^{-1},$ (5.1)

where

${\bf K}=({\bf L}-c{\bf I})({\bf L}+c{\bf
I})^{-1},\quad\mbox{\boldmath${\Upsilon}$}={\bf
I}-2\mbox{\boldmath${\Gamma}$}.$ (5.2)

Since ${\Upsilon}$=$\mbox{\boldmath${\Upsilon}$}^{-1}$ we can rewrite this as

$[-c{\bf I}+\mbox{\boldmath${\Gamma}$}({\bf L}+c{\bf I})]^{-1}=-2({\bf
L}+c{\bf I})^{-1}({\bf I}-\mbox{\boldmath${\Upsilon}$}{\bf
K})^{-1}\mbox{\boldmath${\Upsilon}$}.$ (5.3)

As we are interested in the action of this inverse on ${\bf
h}=\mbox{\boldmath${\Gamma}$}{\bf h}=-\mbox{\boldmath${\Upsilon}$}{\bf h}$ we
can replace the last ${\Upsilon}$ in (5.3) by simply $-{\bf I}$.

Our assumption that the resolvent $(a{\bf I}-{\bf A})^{-1}$ can easily be
computed when $a\notin[0,1]$ allows us to easily compute ${\bf K}$ when
$(1+c)/(1-\sigma)\notin[0,1]$ which is guaranteed to be the case when $c>0$
and $\sigma$ is not on the negative real axis. To bound $\|{\bf K}\|$ we first
consider the special case where $\sigma$ is real with $\sigma>1$ and $c$ is
real and positive with $c\in[1,\sigma]$. Then the norm of ${\bf K}$ as ${\bf
A}$ varies with $0\leq{\bf A}\leq{\bf I}$ is maximized when ${\bf A}=0$ or
${\bf A}={\bf I}$, giving

$\|{\bf
K}\|\leq\max\\{\left|\frac{\sigma-c}{\sigma+c}\right|,\left|\frac{1-c}{1+c}\right|\\},$
(5.4)

and the bound is tightest when $c$ is chosen so the first expression in the
brackets equals the second expression, that is at $c=\sqrt{\sigma}$. With this
value of $c$, and using the fact that $\|\mbox{\boldmath${\Upsilon}$}\|=1$ the
series expansion comes from

$({\bf I}-\mbox{\boldmath${\Upsilon}$}{\bf
K})^{-1}=\sum_{n=0}^{\infty}(\mbox{\boldmath${\Upsilon}$}{\bf K})^{n}\text{
where }\|\mbox{\boldmath${\Upsilon}$}{\bf
K}\|\leq\|\mbox{\boldmath${\Upsilon}$}\|\|{\bf K}\|\leq|u|\text{ with
}u=\frac{\sqrt{\sigma}-1}{\sqrt{\sigma}+1},$ (5.5)

giving

$[-c{\bf I}+\mbox{\boldmath${\Gamma}$}({\bf L}+c{\bf
I})]^{-1}\mbox{\boldmath${\Upsilon}$}=-2({\bf L}+c{\bf
I})^{-1}\sum_{n=0}^{\infty}(\mbox{\boldmath${\Upsilon}$}{\bf K})^{n}.$ (5.6)

So the rate of convergence is $\mu_{3}=|u|$. When $\sigma$ is complex then
with $c$ and $k$ kept at $\sqrt{\sigma}$ and $1/2$ the bound on the norm of
${\bf K}$ needs to be replaced by

$\|{\bf K}\|=\max_{\lambda\in[0,1]}|t(\lambda)|\text{ where
}t(\lambda)=1-\frac{2\sqrt{\sigma}}{\sqrt{\sigma}+1+(\sigma-1)\lambda},$ (5.7)

in which $\lambda$ is a possible eigenvalue of ${\bf A}$. As $t(\lambda)$ is a
fractional linear transformation of $\lambda$ it inscribes a circular arc in
the complex plane joining $t(0)$ with $t(1)=-t(0)$ as $\lambda$ is varied
between 0 and 1, that when extended passes through $t=1$ at $\lambda=\infty$.
The maximum in (5.7) is attained at $\lambda=0$ and $\lambda=1$ if and only if
$t=1$ lies on or outside the circle in the complex $t$ plane centered at $t=0$
and having radius $|t(1)|=-|t(0)|$ , i.e. when $|t(1)|=-|t(0)|\leq 1$. This is
satisfied when $\sigma$ is in the right half of the complex plane and if this
holds then the rate of convergence remains $\mu_{3}=|u|$.

The rate of convergence $\mu_{3}$ is better than $\mu_{1}=|s|$ with
$s=(\sigma-1)/(\sigma+1)$. To see this, we first express $\sqrt{\sigma}$ in
terms of $u$: $\sqrt{\sigma}=(1+u)(1-u)$ and then substitute this back in the
formula for $s$ giving

$s=\frac{2u}{u^{2}+1}.$ (5.8)

So if $u$ and hence $u^{2}$ lie in the unit disk then $|u^{2}+1|\leq 2$
implying that $\mu_{3}\leq\mu_{1}$.

While the rate of convergence $\mu_{3}$ is generally better than the rate of
convergence $\mu_{2}$, i.e. $\mu_{3}<\mu_{2}$ this is not always the case. In
particular, the last series might not converge, i.e. $\mu_{3}=1$, when
$\sigma$ is on the negative real axis, while $\mu_{2}=1$ when $\sigma$ is on a
circle symmetric about the negative real axis intersecting it at the points
$1-1/z^{-}$ and $1-1/z^{+}$ . A comparison of the convergence of the two
series (3.5) and (5.10) for a model example in the theory of composites is
given in [25].

As shown in [15] there is an alternative approach which can be reduced to
essentially the same expansion when ${\bf A}={\bf Q}$. It offers more
flexibility that could result in faster convergence, involving not only a
reference operator ${\bf L}_{0}$, now positive definite but not necessarily
proportional to ${\bf I}$, but also an additional operator ${\bf M}$. Then one
makes a manipulation of the equations into one involving the inverse $({\bf
I}-\mbox{\boldmath${\Upsilon}$}{\bf K})^{-1}$ where now

${\bf K}=[{\bf I}+({\bf L}-{\bf L}_{0}){\bf M}]^{-1}({\bf L}-{\bf
L}_{0}),\quad\mbox{\boldmath${\Upsilon}$}={\bf
M}-\mbox{\boldmath${\Gamma}$}^{\prime}\quad\text{with}\quad\mbox{\boldmath${\Gamma}$}^{\prime}=\mbox{\boldmath${\Gamma}$}(\mbox{\boldmath${\Gamma}$}{\bf
L}_{0}\mbox{\boldmath${\Gamma}$})^{-1}\mbox{\boldmath${\Gamma}$},$ (5.9)

in which the last inverse needs to be taken on the space on which ${\Gamma}$
projects. $\mbox{\boldmath${\Gamma}$}^{\prime}$ can alternatively be defined
through its action: ${\bf E}^{\prime}=\mbox{\boldmath${\Gamma}$}^{\prime}{\bf
P}$ if and only if $\mbox{\boldmath${\Gamma}$}{\bf E}={\bf E}$ and
$\mbox{\boldmath${\Gamma}$}({\bf P}-{\bf L}_{0}{\bf E}^{\prime})=0$. ${\bf M}$
and ${\bf L}_{0}$ need to be such that both ${\bf K}$ and
$\mbox{\boldmath${\Gamma}$}^{\prime}$ can be easily computed. This
manipulation, introduced in [22], and its generalizations are key to
developing alternative formulas for the effective tensors of laminates [14,
32]. Combined with the associated series expansions, it enables the
development of the general theory of exact relations for hierarchical
laminates [8] and for arbitrary microstructures [10] (see also Chapter 17 in
[15] and the book [9]). Exact relations are identities that hold irrespective
of the microstructure provided the material tensors lie on suitable manifolds.

When ${\bf L}_{0}=c{\bf I}$ (now positive definite) and ${\bf M}={\bf I}/2c$
the formulas for ${\bf K}$ and ${\Upsilon}$ , with ${\bf K}$ divided by $2c$
and ${\bf M}$ multiplied by $2c$ match those in (5.2). So we obtain the same
series expansion.

###  5.2 The new series with a faster convergence rate when ${\bf A}$ is a
projection

We now specialize to the case of interest where ${\bf A}={\bf Q}$ is an
orthogonal projection. Our goal is to find a series expansion, with a rate of
convergence at most $\mu_{4}$, which is superior in the sense that $\mu_{4}$
is always less than or equal to $\mu_{2}$ and $\mu_{3}$ (and hence $\mu_{1}$
since $\mu_{2}\leq\mu_{1}$). Then, the resolvent ${\bf R}_{\bf Q}{{\bf
a}}=(a{\bf I}-{\bf Q})^{-1}$ is given by (1.7) and hence ${\bf K}$, given by
(5.2), simplifies to $u(2{\bf Q}-{\bf I})$. Consequently , the series
expansion (5.6) with $c=\sqrt{\sigma}$ becomes

$[{\bf I}-{\mbox{\boldmath${\Gamma}$}}{{\bf Q}}/{z}]^{-1}={\bf
H}\sum_{n=0}^{\infty}u^{n}{\bf D}^{n}\text{ with
}u=\frac{\sqrt{{\sigma}}-1}{\sqrt{{\sigma}}+1},$ (5.10)

where

${\bf H}=2({\bf L}+\sqrt{\sigma}{\bf I})^{-1}=2[{\bf Q}/\sqrt{\sigma}+({\bf
I}-{\bf Q})]/(1+\sqrt{\sigma}),\quad{\bf D}=(2{{\bf Q}}-{\bf I})({\bf
I}-2\mbox{\boldmath${\Gamma}$}),$ (5.11)

and has the same rate of convergence $\mu_{3}=|u|$.

The next step to obtain a series expansion with an improved rate of
convergence is to find one where a fractional linear transformation mapping
the interval $[z^{-},z^{+}]$ in the complex $z$ plane to the negative real
axis naturally enters. The analysis here is based on that given in Chapter 8
of [23]

To start we consider the following linear algebra problem: given
$s_{1},s_{2},s_{3}$ and $t$, solve the matrix equation,

$\begin{pmatrix}J\\\ 0\\\ J_{2}\end{pmatrix}=\begin{pmatrix}E\\\ E_{1}\\\
0\end{pmatrix}-\frac{1}{\underline{z}}\underbrace{\begin{pmatrix}s_{1}^{2}&s_{1}s_{2}&s_{1}s_{3}\\\
s_{1}s_{2}&s_{2}^{2}&s_{2}s_{3}\\\
s_{1}s_{3}&s_{2}s_{3}&s_{3}^{2}\end{pmatrix}}_{{\bf S}}\begin{pmatrix}E\\\
E_{1}\\\ 0\end{pmatrix}-\begin{pmatrix}h\\\ 0\\\ 0\end{pmatrix},$ (5.12)

for $J$ in terms of $E$ and $h$. We will ultimately allow for $s_{1},s_{2}$
and $s_{3}$, that are either real or purely imaginary, chosen with

$s_{3}^{2}=1-s_{1}^{2}-s_{2}^{2},$ (5.13)

to ensure that ${\bf S}$ is a projection matrix, though not self-adjoint in
our application implying that ${\bf S}$ and ${\bf I}-{\bf S}$ project onto
subspaces that are not orthogonal. The significance of (5.12) is that it
corresponds to a problem in the abstract theory of composites, and it enables
us to use the technique of subspace substitution. This technique was
introduced in [15], Section 29.1, including the case where ${\bf S}$ and ${\bf
I}-{\bf S}$ project onto orthogonal subspaces and it was extended in [23],
Sections 7.8 and Chapter 8, to include the case where they project onto non-
orthogonal subspaces. Define ${\cal U}$, ${\cal E}$, and ${\cal J}$ to be the
three subspaces spanned by the three unit vectors

${\bf w}_{0}=\begin{pmatrix}1\\\ 0\\\ 0\end{pmatrix},\quad{\bf
w}_{1}=\begin{pmatrix}0\\\ 1\\\ 0\end{pmatrix},\quad{\bf
w}_{2}=\begin{pmatrix}0\\\ 0\\\ 1\end{pmatrix},$ (5.14)

respectively, so that $\mbox{\boldmath${\Gamma}$}_{i}={\bf w}_{i}\otimes{\bf
w}_{i}$, $i=1,2,3$, are the projections onto ${\cal U}$, ${\cal E}$, and
${\cal J}$ respectively, i.e.

$\mbox{\boldmath${\Gamma}$}_{0}=\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&0\end{pmatrix},\quad\mbox{\boldmath${\Gamma}$}_{1}=\begin{pmatrix}0&0&0\\\
0&1&0\\\
0&0&0\end{pmatrix},\quad\mbox{\boldmath${\Gamma}$}_{2}=\begin{pmatrix}0&0&0\\\
0&0&0\\\ 0&0&1\end{pmatrix}.$ (5.15)

Then (5.12) reduces to

$J{\bf w}_{0}+J_{2}{\bf w}_{2}=\widetilde{{\bf L}}(E{\bf w}_{0}+E_{1}{\bf
w}_{1})-h{\bf w}_{0},\quad\widetilde{{\bf L}}={\bf I}-{\bf S}/\underline{z},$
(5.16)

which is a problem in the abstract theory of composites. More generally, a
problem in the abstract theory of composites takes the form: given ${\bf
E}_{0}\in{\cal U}$, and a source term ${\bf h}$ in ${\cal H}={\cal
U}\oplus{\cal E}\oplus{\cal J}$, and an operator $\widetilde{{\bf L}}$ mapping
${\cal H}$ to ${\cal H}$, find ${\bf J}_{0}\in{\cal U}$, ${\bf E}\in{\cal E}$
and ${\bf J}\in{\cal J}$ such that

${\bf J}_{0}+{\bf J}=\widetilde{{\bf L}}({\bf E}_{0}+{\bf E})-{\bf h}.$ (5.17)

In our case, the subspaces ${\cal U}$, ${\cal E}$, and ${\cal J}$ are clearly
orthogonal, but ${\bf S}$ and ${\bf I}-{\bf S}$ do not generally project onto
orthogonal subspaces when $s_{1},s_{2}$ and $s_{3}$ are not all real.

To find the norm of ${\bf S}$ we consider its action on a possibly complex
vector ${\bf a}$. We have

$|{\bf S}{\bf a}|=|{\bf s}({\bf s}\cdot{\bf a})|\leq|{\bf s}|^{2}|{\bf a}|,$
(5.18)

with equality when ${\bf a}$ is the complex conjugate of ${\bf s}$. Thus ${\bf
S}$ has norm

$\|{\bf S}\|=|s_{1}|^{2}+|s_{2}|^{2}+|s_{3}|^{2},$ (5.19)

and this will surely be greater than or equal to $1$ if (5.13) holds and
$s_{1},s_{2}$ and $s_{3}$ are either real or purely imaginary. For example, if
$s_{1}$ is purely imaginary while $s_{2}$ and $s_{3}$ are purely real then
(5.19) implies

$1=-|s_{1}|^{2}+|s_{2}|^{2}+|s_{3}|^{2}=\|{\bf S}\|^{2}-2|s_{1}|^{2},$ (5.20)

which forces $\|{\bf S}\|$ to be greater than or equal to $1$.

The matrix equation (5.12) is clearly satisfied with
$E_{1}=(s_{1}s_{2})E/(z-s_{2}^{2})$ which gives the so called ”effective
equation”

$J=(1-1/z)E-h,\text{ where }z=\frac{\underline{z}-s_{2}^{2}}{s_{1}^{2}}.$
(5.21)

Solving this last equation for $\underline{z}$ in terms of $z$ gives

$\underline{z}=s_{1}^{2}z+s_{2}^{2}.$ (5.22)

Suppose now that in the extended abstract theory of composites we are
interested in solving the equations

${\bf J}={\bf L}{\bf E}-{\bf
h},\quad\text{with}\quad\mbox{\boldmath${\Gamma}$}{\bf E}={\bf
E},\quad\mbox{\boldmath${\Gamma}$}{\bf J}=0,$ (5.23)

where

${\bf L}={\bf I}-{\bf Q}/z,$ (5.24)

or equivalently in finding the resolvent (1.1) with ${\bf
A}=\mbox{\boldmath${\Gamma}$}{\bf Q}\mbox{\boldmath${\Gamma}$}$. Our
preliminary linear algebra problem shows this is equivalent to solving

$\underbrace{\begin{pmatrix}{\bf J}\\\ 0\\\ {\bf
J}_{2}\end{pmatrix}}_{\underline{{\bf J}}}=\Biggl{[}{\bf
I}-\frac{1}{\underline{z}}\underbrace{\begin{pmatrix}s_{1}^{2}{\bf
Q}&s_{1}s_{2}{\bf Q}&s_{1}s_{3}{\bf Q}\\\ s_{1}s_{2}{\bf Q}&s_{2}^{2}{\bf
Q}&s_{2}s_{3}{\bf Q}\\\ s_{1}s_{3}{\bf Q}&s_{2}s_{3}{\bf Q}&s_{3}^{2}{\bf
Q}\end{pmatrix}}_{\underline{{\bf Q}}}\Biggr{]}\underbrace{\begin{pmatrix}{\bf
E}\\\ {\bf E}_{1}\\\ 0\end{pmatrix}}_{\underline{{\bf
E}}}-\underbrace{\begin{pmatrix}{\bf h}\\\ 0\\\
0\end{pmatrix}}_{\underline{{\bf h}}},$ (5.25)

with $\underline{\mbox{\boldmath${\Gamma}$}}{\bf E}={\bf E}$ and
$\underline{\mbox{\boldmath${\Gamma}$}}{\bf J}=0$, in which

$\underline{\mbox{\boldmath${\Gamma}$}}=\begin{pmatrix}\mbox{\boldmath${\Gamma}$}&0&0\cr
0&{\bf I}&0\cr 0&0&0\end{pmatrix}.$ (5.26)

Specifically, from (5.25) we get

$\begin{pmatrix}{\bf Q}{\bf J}\\\ 0\\\ {\bf Q}{\bf
J}_{2}\end{pmatrix}=\Biggl{[}{\bf
I}-\frac{1}{\underline{z}}\begin{pmatrix}s_{1}^{2}{\bf I}&s_{1}s_{2}{\bf
I}&s_{1}s_{3}{\bf I}\\\ s_{1}s_{2}{\bf I}&s_{2}^{2}{\bf I}&s_{2}s_{3}{\bf
I}\\\ s_{1}s_{3}{\bf I}&s_{2}s_{3}{\bf I}&s_{3}^{2}{\bf
I}\end{pmatrix}\Biggr{]}\begin{pmatrix}{\bf Q}{\bf E}\\\ {\bf Q}{\bf E}_{1}\\\
0\end{pmatrix}-\begin{pmatrix}{\bf Q}{\bf h}\\\ 0\\\ 0\end{pmatrix},$ (5.27)

and

$\begin{pmatrix}({\bf I}-{\bf Q}){\bf J}\\\ 0\\\ ({\bf I}-{\bf Q}){\bf
J}_{2}\end{pmatrix}=\begin{pmatrix}({\bf I}-{\bf Q}){\bf E}\\\ ({\bf I}-{\bf
Q}){\bf E}_{1}\\\ 0\end{pmatrix}-\begin{pmatrix}({\bf I}-{\bf Q}){\bf h}\\\
0\\\ 0\end{pmatrix}.$ (5.28)

Then, from (5.27), (5.12), and (5.21) we obtain

${\bf Q}{\bf J}=(1-1/z){\bf Q}{\bf E}-{\bf Q}{\bf h}=\sigma{\bf Q}{\bf E}-{\bf
Q}{\bf h},$ (5.29)

and (5.28) implies $({\bf I}-{\bf Q}){\bf J}=({\bf I}-{\bf Q}){\bf E}-({\bf
I}-{\bf Q}){\bf h}$. Together they imply (5.23) with ${\bf L}$ given by
(5.24).

We are back at an equivalent problem now involving a new resolvent.
Specifically, since $\underline{{\bf J}}$ and $\underline{{\bf E}}$ lie in
orthogonal spaces, we have

$\underline{{\bf J}}({\bf x})=\underline{{\bf L}}({\bf x})\underline{{\bf
E}}({\bf x})-\underline{{\bf
h}},\quad\underline{\mbox{\boldmath${\Gamma}$}}_{1}\underline{{\bf
E}}=\underline{{\bf
E}},\quad\underline{\mbox{\boldmath${\Gamma}$}}_{1}\underline{{\bf
J}}=0,\quad\text{with}\quad\underline{{\bf L}}={\bf I}-\underline{{\bf
Q}}/\underline{z}.$ (5.30)

Let us now see how this can improve convergence. Applying
$\underline{\mbox{\boldmath${\Gamma}$}}$ to both sides of the first equation
in (5.30) we get

$\underline{{\bf E}}=[{\bf
I}-\underline{\mbox{\boldmath${\Gamma}$}}\underline{{\bf
Q}}/\underline{z}]^{-1}\underline{{\bf h}}=[{\bf
I}-(1-\underline{\sigma})\underline{\mbox{\boldmath${\Gamma}$}}\underline{{\bf
Q}}]^{-1}\underline{{\bf h}},$ (5.31)

in which

$\underline{{\bf Q}}=\begin{pmatrix}s_{1}^{2}{\bf Q}&s_{1}s_{2}{\bf
Q}&s_{1}s_{3}{\bf Q}\\\ s_{1}s_{2}{\bf Q}&s_{2}^{2}{\bf Q}&s_{2}s_{3}{\bf
Q}\\\ s_{1}s_{3}{\bf Q}&s_{2}s_{3}{\bf Q}&s_{3}^{2}{\bf
Q}\end{pmatrix},\quad\underline{\sigma}=1-1/\underline{z}=1+\frac{1}{s_{1}^{2}z+s_{2}^{2}}.$
(5.32)

Note that $\underline{{\bf Q}}$ is a projection operator because both ${\bf
S}$ and ${\bf Q}$ are projections and thus the inverse in (5.31) has exactly
the same form as in (5.10) with $z$ and ${\bf Q}$ being replaced by
$\underline{z}$ and $\underline{{\bf Q}}$ Also
$\underline{\sigma}=(\underline{z}-1)/\underline{z}$ can be re-expressed in
the form

$\underline{\sigma}=\frac{z-z^{-}}{z-z^{+}},$ (5.33)

with

$z^{+}=-s_{2}^{2}/s_{1}^{2},\quad z^{-}=-(1+s_{2}^{2})/s_{1}^{2}.$ (5.34)

Note that $\underline{\sigma}=0$ and $\underline{\sigma}=\infty$ are obtained
by substituting $z=z^{-}$ and $z=z^{+}$ in (5.32). Given real
$1>z^{+}>z^{-}>0$ we need to choose $s_{1}$ and $s_{2}$ so that these
equations are satisfied. This will necessitate complex solutions for $s_{2}$
since otherwise $z^{+}$ will be negative. Explicitly, we have

$s_{1}^{2}=\frac{1}{z^{+}-z^{-}},\quad s_{2}^{2}=\frac{-z^{+}}{z^{+}-z^{-}},$
(5.35)

with $s_{1}$ being real and $s_{2}$ being purely imaginary, and so
$\underline{{\bf Q}}$ is no longer Hermitian implying $\underline{{\bf Q}}$
and ${\bf I}-\underline{{\bf Q}}$ project onto subspaces that are not
orthogonal. This translates to a problem in the extended abstract theory of
composites with a non-orthogonal subspace collection, as introduced in Chapter
8 of [23].

Next, we follow the steps outlined in the previous section, though now the
projection $\underline{{\bf Q}}$ does not have norm 1. We obtain the expansion

$[{\bf I}-{\underline{\mbox{\boldmath${\Gamma}$}}}{\underline{{\bf
Q}}}/\underline{z}]^{-1}=\underline{{\bf
H}}\sum_{n=0}^{\infty}v^{n}\underline{{\bf D}}^{n}$ (5.36)

in powers of

$v=\frac{\sqrt{{\underline{\sigma}}}-1}{\sqrt{{\underline{\sigma}}}+1}=\frac{\sqrt{\frac{z-z^{-}}{z-z^{+}}}-1}{\sqrt{\frac{z-z^{-}}{z-z^{+}}}+1},$
(5.37)

with

$\underline{{\bf H}}=2(\underline{{\bf L}}+\sqrt{\underline{\sigma}}{\bf
I})^{-1}=2\underline{{\bf
Q}/}(\underline{\sigma}+\sqrt{\underline{\sigma}})+2({\bf I}-\underline{{\bf
Q}})/(1+\sqrt{\underline{\sigma}}),\quad\underline{{\bf D}}=(2\underline{{\bf
Q}}-{\bf I})({\bf I}-2\underline{\mbox{\boldmath${\Gamma}$}}),\quad$ (5.38)

Note that the formula for $v$ is a composition of the maps

$v=\frac{w-1}{w+1},\quad
w=\sqrt{\underline{\sigma}},\quad\underline{\sigma}=\frac{z-z^{-}}{z-z^{+}},$
(5.39)

whose action is discussed in Section 3. The desired series expansion is
obtained by substituting the above into the relation

$[{\bf I}-\mbox{\boldmath${\Gamma}$}{\bf Q}/z]^{-1}=\underline{{\bf
G}}^{\dagger}[{\bf I}-\underline{\mbox{\boldmath${\Gamma}$}}\underline{{\bf
Q}}/\underline{z}]^{-1}\underline{{\bf G}},\quad\text{ with }\underline{{\bf
G}}=\begin{pmatrix}{\bf I}\cr 0\cr 0\end{pmatrix},$ (5.40)

giving

$[{\bf I}-\mbox{\boldmath${\Gamma}$}{\bf
Q}/z]^{-1}=\frac{2}{\sqrt{\underline{\sigma}}+1}\sum_{n=0}^{\infty}v^{n}\underline{{\bf
G}}^{\dagger}\underline{{\bf H}}\underline{{\bf D}}^{n}\underline{{\bf G}}.$
(5.41)

As this is the expansion of a self-adjoint operator in powers of $v$ we deduce
that the terms
$\underline{{\bf G}}^{\dagger}\underline{{\bf H}}\underline{{\bf
D}}^{n}\underline{{\bf G}}$ must be self-adjoint even though $\underline{{\bf
D}}$ is not. Since $[{\bf I}-\mbox{\boldmath${\Gamma}$}{\bf Q}/z]^{-1}$ has
potentially singularities anywhere in the interval $[z^{-},z^{+}]$,
corresponding to $v$ being on the unit circle we deduce that the radius of
convergence is at least $1$ so the rate of convergence is at least
$\mu_{4}=|v|$. To compare the rates of convergence bounds $\mu_{2}$ and
$\mu_{4}$ it is helpful to first express both in terms of
$w=\sqrt{\underline{\sigma}}$:

$\mu_{2}=\left|\frac{w^{2}-1}{w^{2}+1}\right|,\quad\mu_{4}=\left|\frac{w-1}{w+1}\right|.$
(5.42)

Thus the convergent rate of the first series with a given value
$\underline{\sigma}=w$ will be the same as the second when $\sigma=w^{2}$. A
similar analysis as that which lead to (5.8) now implies

$\mu_{2}=\frac{2|v|}{|v^{2}+1|},\quad\mu_{4}=|v|.$ (5.43)

Clearly if $v$ is in the unit disk then so is $v^{2}$ and $|v^{2}+1|$ is at
most 2, equaling 2 when $v=\pm 1$. As the unit disk is the image of the
complex $z$-plane, excluding the cut $z\in[z^{-},z^{+}]$ we conclude that
$\mu_{4}\leq\mu_{2}$. When $v$ is small then $\mu_{4}\approx 2\mu_{4}$, both
being very large. We do not have a proof that $\mu_{4}\leq\mu_{3}$, although
we suspect that to be the case. The convergence has been tested by Moulinec
and Suquet (private communication) in a model example of a conducting
composite and the convergence based on the series (5.38) is better than that
of the series (5.10): see figure 8.8 in [23].

Naturally, computing the action of $\underline{{\bf D}}$ on a field is more
involved than computing the action of ${\bf D}$ on a field. So, if $z^{-}$ is
close to zero and $z^{+}$ is close to 1, it is probably best to work with the
series (5.10) rather than (5.38).

###  5.3 What if $z^{+}$ and $z^{+}$ are only known approximately?

We may not know $z^{+}$ and $z^{+}$ exactly but rather just have estimates
$z^{+}_{e}$ and $z^{+}_{e}$ for them. The treatment here follows that
discussed in Section 8.8 of [23]. Using the estimates gives

$v=\left(\frac{\sqrt{\frac{z-z_{e}^{-}}{z-z_{e}^{+}}}-1}{\sqrt{\frac{z-z_{e}^{-}}{z-z_{e}^{+}}}+1}.\right)$
(5.44)

The interesting case is when $z^{+}>z^{+}_{e}$ or $z^{-}_{e}>z^{-}$ since
otherwise $z^{-}_{e}{\bf I}\leq\mbox{\boldmath${\Gamma}$}{\bf
A}\mbox{\boldmath${\Gamma}$}\leq z^{+}_{e}{\bf I}$ and all the previous
results hold with $z^{+}$ and $z^{-}$ replaced with $z^{+}_{e}$ and
$z^{-}_{e}$. Just looking at the case where $z^{+}>z^{+}_{e}$ and
$z^{-}_{e}>z^{-}$ the points $z^{-}$ and $z^{+}$ map to

$v^{-}=\left(\frac{\sqrt{\frac{z^{-}-z_{e}^{-}}{z^{-}-z_{e}^{+}}}-1}{\sqrt{\frac{z^{-}-z_{e}^{-}}{z^{-}-z_{e}^{+}}}+1}\right),\quad
v^{+}=\left(\frac{\sqrt{\frac{z^{+}-z_{e}^{-}}{z^{+}-z_{e}^{+}}}-1}{\sqrt{\frac{z^{+}-z_{e}^{-}}{z^{+}-z_{e}^{+}}}+1}\right).$
(5.45)

To the leading order in $z^{-}-z_{e}^{-}$ and $z^{+}-z_{e}^{+}$ approximate to

$v^{-}\approx\sqrt{\frac{z_{e}^{-}-z^{-}}{z^{+}-z^{-}}}-1,\quad v^{+}\approx
1-\sqrt{\frac{z^{+}-z_{e}^{+}}{z^{+}-z^{-}}}.$ (5.46)

Assuming that the spectrum of ${\bf A}$ includes both $z^{-}$ and $z^{+}$, the
rate of convergence will be the greater of $|v/v^{+}|$ and $|v/v^{-}|$. Since
(5.46) implies $|v^{-}|\approx 1$ when $z^{-}-z_{e}^{-}$ is small and
$|v^{-}+|\approx 1$ when $z^{+}-z_{e}^{+}$ is small the rate of convergence
will be close to that obtained when we know $z^{+}$ and $z^{-}$ exactly.

A natural way of obtain these estimates $z_{e}^{-}$ and $z_{e}^{+}$ is via the
Rayleigh Ritz procedure. Thus, we pick appropriate finite dimensional
subspaces ${\cal S}^{-}\subset{\mbox{\boldmath${\cal E}$}}$ and ${\cal
S}^{+}\subset{\mbox{\boldmath${\cal E}$}}$ and take

$\displaystyle z^{-}_{e}$ $\displaystyle=$ $\displaystyle\min_{{\bf a}\in{\cal
S}^{-}\,\,|{\bf a}|=1}({\bf P}^{\dagger}{\bf A}{\bf P}{\bf a},{\bf a}),$
$\displaystyle z^{+}_{e}$ $\displaystyle=$ $\displaystyle\max_{{\bf a}\in{\cal
S}^{+},\,\,|{\bf a}|=1}({\bf P}^{\dagger}{\bf A}{\bf P}{\bf a},{\bf a}),$
(5.47)

where $(\cdot,\cdot)$ denotes the inner product in ${\cal E}$ which we take to
be the inner product in our Hilbert or vector space ${\cal H}$. If, for
example, the spectrum was discrete then ideally there should be a field ${\bf
a}\in{\cal S}^{-}$ that is a good approximation to some field in the
eigenspace of $\mbox{\boldmath${\Gamma}$}{\bf A}\mbox{\boldmath${\Gamma}$}$
having the lowest eigenvalue, and ideally there should be a field ${\bf
a}\in{\cal S}^{+}$ that is a good approximation to some field in the
eigenspace of $\mbox{\boldmath${\Gamma}$}{\bf A}\mbox{\boldmath${\Gamma}$}$
having the maximum eigenvalue.

##  6 Conclusion

We have presented a novel approach to calculating a certain class of
resolvents via series expansions. These series expansions provide a generally
rapidly converging way of computing the action of these resolvents on fields
via an iterative method. While the method assumes the spectrum of ${\bf
P}^{\dagger}{\bf Q}{\bf P}$ lies inside a known, or at least approximately
known, interval $[z^{-},z^{+}]$, it seems likely that adaptive methods could
be developed that do not require knowledge of $z^{-}$ and $z^{+}$. The
performance of the method improves upon some methods, and likely upon all the
alternate methods presented here, at least in the asymptotic limit where one
is close to convergence. However, it remains to provide bounds on the error
incurred when one truncates the series expansion. Clearly, explicit numerical
investigations, beyond those of Moulinec and Suquet (private communication,
2015) as summarized in figure 8.8 of [23] need to be made. These should
compare the usefulness of the method with other available methods, both in
terms of speed of convergence and computer memory required. In particular, it
remains to compare the method presented here with more established methods,
such as the conjugate gradient method (with $z$ real) or the biconjugate
gradient method, biconjugate gradient stabilized method, or conjugate gradient
squared method (with $z$ complex). When $z$ is real the condition number of
the operator $(z{\bf I}-{\bf P}^{\dagger}{\bf Q}{\bf P})$ we are inverting is
$\kappa=\underline{\sigma}=(z-z^{-})/(z-z^{+})$ when $z>z^{+}$ and
$\kappa=1/\underline{\sigma}$ when $z<z^{-}$. So the convergence rate of our
new method

$|v|=\left|\frac{\sqrt{\underline{\sigma}}-1}{\sqrt{\underline{\sigma}}+1}\right|=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}$
(6.48)

matches that of the conjugate gradient method. Comparisons of convergence
between the conjugate gradient method and the basic scheme [26] we discussed
in section 3 have been made for conducting composites [31] and elastic
composites [2]. As far as memory goes the new method is likely to be more
demanding than most other methods because $\underline{{\bf Q}}$ and
$\underline{\mbox{\boldmath${\Gamma}$}}$ act in spaces that have three times
the dimensionality of ${\cal H}$ when ${\cal H}$ is finite dimensional.

Even if the conjugate gradient, biconjugate gradient, biconjugate gradient
stabilized, and conjugate gradient squared methods perform better than the
method ultimately presented here, the idea of subspace substitution with non-
orthogonal subspaces may have application in other areas of mathematics and
numerical analysis in particular.

## Acknowledgments

Matti Schneider is thanked for many helpful comments on the manuscript and
Herve Moulinec and Pierre Suquet are thanked for additional comments. The
author is grateful to the National Science Foundation for support through
Research Grant DMS-2107926.

## References

  * [1] A. Braides, I. Fonseca, and G. Leoni, A-quasiconvexity: Relaxation and homogenization, ESIAM: COCV, 5 (2000), pp. 539–577.
  * [2] S. Brisard and L. Dormieux, FFT-based methods for the mechanics of composites: A general variational framework, Computational Materials Science, 49 (2010), pp. 663–671.
  * [3] W. F. Brown, Jr., Solid mixture permittivities, Journal of Chemical Physics, 23 (1955), pp. 1514–1517.
  * [4] B. Dacorogna, Direct Methods in the Calculus of Variations, vol. 78 of Applied Mathematical Sciences, , Springer-Verlag, Berlin, Germany / Heidelberg, Germany / London, UK / etc., second ed., 2007.  Series editors: S. S. Antman, J. E. Marsden, and L. Sirovich.
  * [5] J. Douglas, Jr and H. Rachford, Jr, On the numerical solution of heat conduction problems in two and three space variables, Transactions of the American Mathematical Society, 82 (1956), pp. 421–439.
  * [6] D. J. Eyre and G. W. Milton, A fast numerical scheme for computing the response of composites using grid refinement, European Physical Journal. Applied Physics, 6 (1999), pp. 41–47.
  * [7] P. Giselsson and S. Boyd, Linear convergence and metric selection for Douglas-Rachford splitting and admm, IEEE Transactions on Automatic Control, 62 (2017), pp. 532–544.
  * [8] Y. Grabovsky, Exact relations for effective tensors of polycrystals. I: Necessary conditions, Archive for Rational Mechanics and Analysis, 143 (1998), pp. 309–329.
  * [9] , Composite Materials: Mathematical Theory and Exact Relations, IOP Publishing, Bristol, UK, 2016.
  * [10] Y. Grabovsky, G. W. Milton, and D. S. Sage, Exact relations for effective tensors of composites: Necessary conditions and sufficient conditions, Communications on Pure and Applied Mathematics (New York), 53 (2000), pp. 300–353.
  * [11] E. Kröner, Statistical Continuum Mechanics: Course Held at the Department of General Mechanics, Udine, October, 1971, vol. 92 of International Center for Mechanical Sciences, Courses and Lectures, Springer-Verlag, Berlin / Heidelberg / London / etc., 1972, pp. 109–120.
  * [12] J. C. Michel, H. Moulinec, and P. M. Suquet, A computational method based on augmented Lagrangians and Fast Fourier Transforms for composites with high contrast, Computer Modeling in Engineering and Sciences, 1 (2000), pp. 79–88.
  * [13] , A computational scheme for linear and nonlinear composites with arbitrary phase contrast, International Journal for Numerical Methods in Engineering, 52 (2001), pp. 139–160.
  * [14] G. W. Milton, On characterizing the set of possible effective tensors of composites: The variational method and the translation method, Communications on Pure and Applied Mathematics (New York), 43 (1990), pp. 63–125.
  * [15] , The Theory of Composites, vol. 6 of Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, Cambridge, UK, 2002.  Series editors: P. G. Ciarlet, A. Iserles, Robert V. Kohn, and M. H. Wright.
  * [16] , A unifying perspective on linear continuum equations prevalent in physics. Part I: Canonical forms for static and quasistatic equations.  Available as arXiv:2006.02215 [math.AP]., 2020.
  * [17] , A unifying perspective on linear continuum equations prevalent in physics. Part II: Canonical forms for time-harmonic equations. Available as arXiv:2006.02433 [math-ph]., 2020.
  * [18] , A unifying perspective on linear continuum equations prevalent in physics. Part III: Canonical forms for dynamic equations with moduli that may, or may not, vary with time. Available as arXiv:2006.02432 [math-ph], 2020.
  * [19] , A unifying perspective on linear continuum equations prevalent in physics. Part IV: Canonical forms for equations involving higher order gradients. Available as arXiv:2006.03161 [math-ph]., 2020.
  * [20] , A unifying perspective on linear continuum equations prevalent in science. Part V: resolvents; bounds on their spectrum; and their Stieltjes integral representations when the operator is not selfadjoint. Available as arXiv:2006.03162 [math-ph], 2020.
  * [21] , A unifying perspective on linear continuum equations prevalent in science. Part VI: rapidly converging series expansions for their solution.  Available as arXiv:2006.11388 [math-ph], 2020.
  * [22] G. W. Milton and K. M. Golden, Representations for the conductivity functions of multicomponent composites, Communications on Pure and Applied Mathematics (New York), 43 (1990), pp. 647–671.
  * [23] G. W. Milton (editor), Extending the Theory of Composites to Other Areas of Science, Milton–Patton Publishers, P.O. Box 581077, Salt Lake City, UT 85148, USA, 2016.  Available from bookbaby.com.
  * [24] H. Moulinec and F. Silva, Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials, International Journal for Numerical Methods in Engineering, 97 (2014), pp. 960–985.
  * [25] H. Moulinec, P. Suquet, and G. W. Milton, Convergence of iterative methods based on Neumann series for composite materials: theory and practice, International Journal for Numerical Methods in Engineering, 114 (2018), pp. 1103–1130.
  * [26] H. Moulinec and P. M. Suquet, A fast numerical method for computing the linear and non-linear properties of composites, Comptes rendus des Séances de l’Académie des sciences. Série II, Tome 318 (1994), pp. 1417–1423.
  * [27] M. Schneider, On the Barzilai-Borwein basic scheme in FFT-based computational homogenization, International Journal for Numerical Methods in Engineering, 118 (2019), pp. 482–494.
  * [28] , A review of nonlinear FFT-based computational homogenization methods, Acta Mechanica, 232 (2021), pp. 2051–2100.
  * [29] M. Schneider, D. Wicht, and T. Böhlke, On polarization-based schemes for the FFT-based computational homogenization of inelastic materials, Computational Mechanics, 64 (2021), pp. 1073–1095.
  * [30] G. Strang, Introduction to Applied Mathematics, Wellesley–Cambridge Press, Wellesley, Massachusetts, 1986.
  * [31] J. Zeman, J. Vondřejc, J. Nová, and I. Marek, Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients, Journal of Computational Physics, 229 (2010), pp. 8065–8071.
  * [32] V. V. Zhikov, Estimates for the homogenized matrix and the homogenized tensor, Uspekhi Matematicheskikh Nauk = Russian Mathematical Surveys, 46 (1991), pp. 49–109.  English translation in Russ. Math. Surv. 46(3):65–136 (1991).