An extremal problem arising in the dynamics of two-phase materials that directly reveals information about the internal geometry

Many systems have responses that are nonlocal in time as this is naturally a consequence of the fact that it takes time for subelements of the system to respond. Usually, this leads one to either: (1) examine the response at each, or one or more, frequencies as the convolution in time characterizing the response of the system becomes a simple multiplication in the frequency domain; or (2) examine the response to a delta function or Heaviside function as this directly reveals the integral kernel characterizing the response. In this paper we show that desired information about the system can be directly obtained from selectively designed input signals that are neither at constant frequency, nor delta or Heaviside functions.

Our initial motivation comes from the work [31, 33] where we derived microstructure-independent bounds on the viscoelastic response at a given time $t$ of two-phase periodic composites (in antiplane shear) with prescribed volume fractions $f_{1}$ and $f_{2}=1-f_{1}$ of the phases and with an applied average stress or strain prescribed as a function of time. We found that the bounds were sometimes extremely tight at particular times $t=t_{0}$: see Figure 1. This was quite a surprise because the response of each phase is nonlocal in time, yet somehow this response is untangled at these particular times. Thus, the bounds could be used in an inverse fashion to determine the volume fractions from measurements at time $t_{0}$. Determining volume fractions of phases is important in the oil industry, where one wants to know the proportions occupied by oil and water in the rock, to detecting breast cancer, to assessing the porosity of sea-ice and other materials, and even to determining the volume of holes in swiss cheese. While our bounds were very tight at specific times in some examples, they were far from tight at all times in other examples: see Figure 2.

At that time it was totally unclear as to whether appropriate input signals could produce the desired tight bounds at a specific time and, if so, what algorithm should be followed to design these input signals. The primary goal of this paper is to address this problem in the case where the input function has a finite number of frequencies. In particular, for the example of Figure 2, our algorithm produces the much tighter bounds of Figure 3. It is an open question as to whether one can find smooth input signals, containing a continuum of frequencies, such that the response $\mathop{\rm Re}\nolimits[v(t_{0})]$ of the material at a specific moment of time $t_{0}$ is totally measure independent, while $\mathop{\rm Re}\nolimits[v(t)]$ has a smooth dependence on $t$, with $\mathop{\rm Re}\nolimits[v(t)]\to 0$ when $t\to-\infty$. The example of Figure 1 suggests that it may be possible.

We emphasize that our results are applicable not just to determining the volume fractions of the phases in a two-phase composite but also determining the volume and shape of an inclusion in a body from exterior boundary measurements. This is shown in Sections 10 and 11. It is a classical and important inverse problem with a long history and many contributions: see [10, 26, 27, 47, 46, 22, 1, 19, 9, 11, 4, 3, 25, 44, 28] and references therein.

A secondary goal of this paper is to solve an accompanying mathematical approximation problem which we now outline, and which is essential to achieve the primary goal.

Here we formulate a problem, not just relevant to systems with a nonlocal response, but potentially to other application where Markov functions arise. Such functions are Cauchy transforms of positive measures with compact support on the real line. Some authors may suitably call them Herglotz functions, or Nevanlinna functions, or Stieltjes transforms.

Suppose $F_{\mu}(z)$ is a Markov function having the integral representation

$$F_{\mu}(z)=\int_{-1}^{1}\frac{d\mu(\lambda)}{\lambda-z},$$

(2.1)

where the Borel measure $\mu$ is positive with unit mass:

$$\int_{-1}^{1}d\mu(\lambda)=1.$$

(2.2)

Given $m$ (possibly complex) points $z_{1},z_{2},\ldots,z_{m}$ not belonging to the interval $[-1,1]$, we are interested in finding complex constants $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$ such that

$$\sum_{k=1}^{m}\alpha_{k}F_{\mu}(z_{k})\approx 1$$

(2.3)

for all probability measures $\mu$. Optimal bounds correlating the possible values of the $m$-tuple $(F_{\mu}(z_{1}),F_{\mu}(z_{2}),\ldots,F_{\mu}(z_{m}))$ as $\mu$ varies over all probability measures are well known, as derived from the well charted analysis of the Nevanlinna-Pick interpolation problem [29]. Indeed, the nonlinear constraints among the values $F_{\mu}(z_{1}),F_{\mu}(z_{2}),\ldots,F_{\mu}(z_{m})$, and standard convexity theory provide optimal bounds on the range of the left hand side of (2.3) for given constants $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$, see [29] for details. But this is not our main concern.

We would rather like to choose $m$ points $z_{1},z_{2},\ldots,z_{m}$, and find associated constants $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$, for every prescribed integer $m$, having the property

$$\sup_{\mu}\left|\sum_{k=1}^{m}\alpha_{k}F_{\mu}(z_{k})-1\right|\leq\epsilon_{m}$$

(2.4)

for some bound $\epsilon_{m}$ subject to $\epsilon_{m}\to 0$ as $m\to\infty$. The geometry of the locus of these points is obviously essential and it will be detailed in the sequel. The faster the convergence, the better.

Since we deal with probability measures, condition (2.4) is equivalent to

$$\sup_{\lambda\in[-1,1]}\left|\sum_{k=1}^{m}\frac{\alpha_{k}}{\lambda-z_{k}}-1\right|\leq\epsilon_{m}.$$

(2.5)

And this is good news. We seek the minimal deviation from zero on the interval $[-1,1]$ of a rational function $R(z)$ satisfying $R(\infty)=1$ and possessing simple poles at the points $z_{1},\ldots,z_{m}$. Or equivalently, denoting $q(z)=(z-z_{1})(z-z_{2})\ldots(z-z_{m})$ and $w(\lambda)=|q(\lambda)|^{-1}$ we aim at finding the minimal deviation from zero of a monic polynomial $p$ of degree $m$, with respect to the weighted norm $\|p\ w\|_{\infty}=\sup_{\lambda\in[-1,1]}|p(\lambda)w(\lambda)|.$

Both perspectives align to well-known classical studies in approximation theory. The first one is an extremal problem in rational approximation with prescribed poles, a subject going back at least to Walsh [51]. A great deal of information in this respect was systematized in Walsh’s book [52]. The second approach is a genuine weighted Chebyshev approximation problem, and here we are on solid ground. First, note that the functions

$$w(\lambda),\lambda w(\lambda),\lambda^{2}w(\lambda),\ldots$$

(2.6)

form a Chebyshev system on the interval $[-1,1]$, that is, they are linearly independent and any linear combination of $w(\lambda),\lambda w(\lambda),\lambda^{2}w(\lambda),\ldots,\lambda^{m}w(\lambda)$ has at most $m$ zeros in $[-1,1]$. Even more, a stronger so-called Markov property of this system of functions holds. The classical Chebyshev approximation in the uniform norm theorem has an analog for such non-orthogonal bases [23, 29]. To be more precise, there exists a unique monic polynomial $p$ of degree $m$ minimizing the norm $\|p\ w\|_{\infty}$: this polynomial is characterized by the fact that $|p(\lambda)|$ attains its maximal value at $m+1$ points, and the sign of $p(\lambda)$ alternates there, see also [34]. In case $w(\lambda)=1$, the optimal polynomial is of course the normalized Chebyshev polynomial of the first kind: $p(\lambda)=\frac{T_{m}(\lambda)}{2^{m}},\ T_{m}(\cos x)=\cos(mx),\,m\geq 0.$ The constructive aspects of weighted Chebyshev approximation are rather involved, see for instance the early works of Werner [54, 55, 56]. In the same vein, the asymptotics of the optimal bound of our minimization problem inherently involves potential theory or operator theory concepts. We cite for a comparison basis a few remarkable results of the same flavor [16, 45, 6].

Without seeking sharp bounds and guided by the specific applications we aim at, we propose a compromise and relaxation of our extremal problem:

$$\inf_{p}\|p\ w\|_{\infty}\leq\|w\|_{\infty}\inf_{p}\|p\|_{\infty}.$$

(2.7)

At this point we can invoke Chebyshev original theorem and his polynomial $T_{m}$, obtaining in this way the benefit, very useful for applications, of computing in closed form the residues $\alpha_{k}$. Details and some ramifications will be given in the Section 4.

Without going into the specific details, as these will be provided later in Section 7, 8 and 9, in many linear systems with an input function $u(t)$ varying with time $t$, of the form

$$u(t)=\sum_{k=1}^{m}\beta_{k}e^{-i\omega_{k}(t-t_{0})},$$

(3.1)

where the $\omega_{k}$ are a set of (possibly complex) frequencies, and $t_{0}$ is a given time, the output function $v(t)$ takes the form

$$v(t)=\sum_{k=1}^{m}\alpha_{k}a_{0}F_{\mu}(z(\omega_{k}))e^{-i\omega_{k}(t-t_{0})},$$

(3.2)

in which the function $F_{\mu}(z)$ is given by (2.1),

$$\alpha_{k}=\beta_{k}c(\omega_{k}),$$

(3.3)

and the functions $z(\omega)$ and $c(\omega)$ depend on $\omega$ in some known way: $z=z(\omega)$ and $c=c(\omega)$. The real constant $a_{0}>0$ and the unknown measure $d\mu$ depend on the system. In our viscoelasticity study [33] the connection with Markov functions comes from the fact that the effective shear modulus $G_{*}(\omega)$, that relates the average stress to the average strain at frequency $\omega$, as a function of the shear moduli $G_{1}(\omega)$ and $G_{2}(\omega)$ of the two phases, has the property that $[(G_{*}/G_{1})-1]/(2f_{1})$, in which $f_{1}$ is the volume fraction of phase 1, is a Markov function of $z=(G_{1}+G_{2})/(G_{2}-G_{1})$ taking the form (2.2) [7, 35, 15].

Henceforth, we adopt the notational simplification

$$f(z)=F_{\mu}(z).$$

Thus, at time $t=t_{0}$, the output function is

$$v(t_{0})=a_{0}\sum_{k=1}^{m}\alpha_{k}f(z_{k})\quad\text{with}\quad z_{k}=z(\omega_{k}),$$

(3.4)

and we seek an input signal so that the output $v(t_{0})$ is almost system independent with $v(t_{0})\approx a_{0}$. So, by measuring $v(t_{0})$ we can determine the system parameter $a_{0}$. In the viscoelastic problem that we studied [31, 33], $a_{0}$ is the volume fraction $f_{1}$ (see also [7]) and it is useful to be able to determine this from indirect measurements. Typically, one may assume the frequencies $\omega_{k}$ have a positive imaginary part so that the input signal $u(t)$ is essentially zero in the distant past. In (3.1) one could just take a signal with $m-1$ frequencies $\omega_{k}$, $k=1,2,\ldots,m-1$. Then, we have

$$v(t_{0})+\alpha_{m}a_{0}f(z(\omega_{m}))=a_{0}\sum_{k=1}^{m}\alpha_{k}f(z_{k})\approx a_{0}\quad\text{with}\quad z_{k}=z(\omega_{k}).$$

(3.5)

Then, if $a_{0}$ is known, a measurement of $v(t_{0})$ will allow us to estimate the output $a_{0}f(z(\omega_{m}))e^{-i\omega_{m}(t-t_{0})}$ at a desired (possibly real) frequency $\omega_{m}$ given the time harmonic input $e^{-i\omega_{m}(t-t_{0})}$.

It is often the situation, such as in the viscoelastic problem, that only the real part of $v(t)$ has a direct physical significance and, hence, one might want to find constants $\alpha_{k}$ such that, say,

	$\displaystyle 2\mathop{\rm Re}\nolimits[v(t_{0})]$	$\displaystyle=a_{0}\left(\sum_{k=1}^{m}\alpha_{k}f(z_{k})+\sum_{k=1}^{m}\overline{\alpha_{k}}\overline{f(z_{k})}\right)$
		$\displaystyle=a_{0}\left(\sum_{k=1}^{m}\alpha_{k}f(z_{k})+\sum_{k=1}^{m}\overline{\alpha_{k}}f(\overline{z_{k}})\right)\approx a_{0},$		(3.6)

where the overline denotes complex conjugation. This, again, reduces to a problem of the form (2.3) where, after renumbering, the complex values of $z_{k}$ come in pairs, $z_{k}$ and $z_{k+1}=\overline{z_{k}}$ and we may take $\alpha_{k+1}=\overline{\alpha_{k}}$ so that the left hand side of (2.3) is real.

It may be the case that the first $n$ moments of the probability measure $d\mu$ are known,

$$M_{\ell}=\int_{-1}^{1}\lambda^{\ell}d\mu(\lambda),\quad\ell=1,2,\ldots,n,$$

(3.7)

in addition to the zeroth moment $M_{0}=1$ and that $m$ (possibly complex) points $z_{1},z_{2},\ldots,z_{m}$ not on the interval $[-1,1]$ are given. We then may seek complex constants $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$ and $\gamma_{0},\gamma_{1},\gamma_{2},\ldots\gamma_{n}$, with say $\gamma_{n}=1$, such that

$$\sum_{k=1}^{m}\alpha_{k}f(z_{k})\approx\sum_{\ell=0}^{n}\gamma_{\ell}M_{\ell}$$

(3.8)

for all probability measures $\mu$ with the prescribed moments. We will treat this problem in Section 3. We can use these results to determine an approximate linear relation among the $n$ moments if $v(t_{0})$ is measured. This may be used to estimate one moment if the rest are known. This can be useful when the moments have an important physical significance: in the viscoelastic problem, for instance, $M_{1}$ depends only on the volume fraction $f_{1}$ if one assumes that the composite has sufficient symmetry to ensure that its response remains invariant as the material is rotated [7]. So, incorporating the moment $M_{1}$ and measuring the response at time $t_{0}$, then, allows us to obtain tighter bounds on $f_{1}$ as shown in [31, 33].

Mutatis mutandis, we may seek complex constants $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$ and $\gamma_{1},\gamma_{2},\ldots\gamma_{n}$ (each constant depending both on $m$ and $n$) such that

$$\inf_{\bf A}\left\|\sum_{k=1}^{m}\alpha_{k}[{\bf A}-z_{k}{\bf I}]^{-1}-\sum_{\ell=0}^{n}\gamma_{\ell}{\bf A}^{\ell}\right\|\leq\epsilon_{m}^{(n)},$$

(3.9)

where the infimum is over all self adjoint operators ${\bf A}$ with spectrum in $[-1,1]$ and for fixed $n$, $\epsilon_{m}^{(n)}\to 0$ as $m\to\infty$. Our analysis extends easily to treat this problem too. The relevance is that in many linear systems with an input field ${\bf u}(t)$ varying with time $t$, of the form

$${\bf u}(t)=\sum_{k=1}^{m}\beta_{k}e^{-i\omega_{k}(t-t_{0})}{\bf u}_{0},$$

(3.10)

the output field ${\bf v}(t)$ takes the form

$${\bf v}(t)=\sum_{k=1}^{m}\alpha_{k}e^{-i\omega_{k}(t-t_{0})}a_{0}[{\bf A}-z(\omega_{k}){\bf I}]^{-1}{\bf u}_{0}\quad\text{with}\quad\alpha_{k}=\beta_{k}c(\omega_{k}),$$

(3.11)

where the real constant $a_{0}$ and the self-adjoint operator ${\bf A}$ characterize the response of the system, and the system parameters $z(\omega)$ and $c(\omega)$ depend on the frequency $\omega$ in some known way. Then, the bound (3.9) implies

$$\left|{\bf v}(t_{0})-a_{0}\sum_{\ell=0}^{n}\gamma_{\ell}{\bf A}^{\ell}{\bf u}_{0}\right|\leq a_{0}\epsilon_{m}^{(n)}\left|{\bf u}_{0}\right|.$$

(3.12)

The present section contains the main result which provides the theoretical foundation of our explorations. As explained in the introduction, we try to balance the computational accessibility and simplicity with the loss of sharp bounds. A few comments about the versatility of the following theorem are elaborated after its proof.

Now we may write

$$\sum_{k=1}^{m}\frac{\alpha_{k}}{\lambda-z_{k}}=\frac{p(\lambda)}{q(\lambda)}=R(\lambda),$$

(4.8)

where $q(\lambda)$ is the known monic polynomial

$$q(\lambda)=\prod_{j=1}^{m}(\lambda-z_{j})$$

(4.9)

of degree $m$, and $p(\lambda)$ is a polynomial of degree at most $m-1$ that remains to be determined. The constants $\alpha_{k}$ can then be identified with the residues at the poles $\lambda=z_{k}$ of $R(\lambda)$:

$$\alpha_{k}=\frac{p(z_{k})}{\prod_{j\neq k}(z_{k}-z_{j})}.$$

(4.10)

Then, the problem becomes one of choosing $p(\lambda)$ such that

$$\sup_{\lambda\in[-1,1]}\left|\frac{p(\lambda)}{q(\lambda)}-1\right|=\sup_{\lambda\in[-1,1]}\frac{\left|p(\lambda)-q(\lambda)\right|}{|q(\lambda)|}$$

(4.11)

is close to zero. Clearly, the problem is now one of polynomial approximation of the monic polynomial $q(\lambda)$ of degree $m$ by the polynomial $p(\lambda)$. A natural choice is to take

$$p(\lambda)=q(\lambda)-T_{m}(\lambda)/2^{m-1},$$

(4.12)

where $T_{m}(\lambda)/2^{m-1}$ is the Chebyshev polynomial $T_{m}(\lambda)$ of degree $m$, normalized to be monic. This choice minimizes the sup-norm of $|p(\lambda)-q(\lambda)|$ over the interval $\lambda\in[-1,1]$ and

$$|p(\lambda)-q(\lambda)|=|T_{m}(\lambda)/2^{m-1}|\leq 1/2^{m-1}$$

(4.13)

provides a bound on the numerator in (4.11). To bound the denominator, we have

$$|q(\lambda)|=\prod_{k=1}^{m}|\lambda-z_{k}|\geq\prod_{k=1}^{m}d(z_{k}),$$

(4.14)

where $d(z_{k})$ is given by (4.1). Using (4.2) and the bounds (4.13) and (4.14) we see that (2.4) is satisfied with $\epsilon^{m}=2/(2d_{\rm min})^{m}$. Finally, with $p(\lambda)$ given by (4.12) we see that the residues $\alpha_{k}$ at the poles $\lambda=z_{k}$ of $g(\lambda)$, given by (4.10) correspond to those given by (4.3).

By taking the natural logarithm, we are led to enforce the condition

$$\limsup_{n}\sup_{\lambda\in[-1,1]}\frac{1}{n}\sum_{j=1}^{n}\ln\frac{1}{|\lambda-z_{j}(n)|}<\ln 2.$$

In other terms, an evenly distributed probability mass on the points $z_{1}(n),\ldots,z_{n}(n)$ should have its logarithmic potential asymptotically bounded from above by a prescribed constant, on the interval $[-1,1]$. Again, this turns out to be a rather typical problem of approximation theory, at least when restricting the poles of $q_{n}$ to belong to some Jordan curve surrounding $[-1,1]$. A natural choice being an ellipse with foci at $\pm 1$, see also [45, 6].

Here we assume that the first $n$ moments $M_{1},M_{2},\ldots,M_{n}$ of the probability measure $d\mu$, given by (2.1), are known, in addition to $M_{0}=1$ and that $m$ (possibly complex) points $z_{1},z_{2},\ldots,z_{m}$ not on the interval $[-1,1]$ are given. We seek complex constants $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$ and $\gamma_{1},\gamma_{2},\ldots\gamma_{n}$, with say $\gamma_{n}=1$ such that

$$\left|\sum_{k=1}^{m}\alpha_{k}f(z_{k})-\sum_{\ell=0}^{n}\gamma_{\ell}M_{\ell}\right|$$

(5.1)

is small for all probability measures $\mu$ with the prescribed $n$ moments. The analysis proceeds as before, only now we introduce the polynomial

$$r(\lambda)=\sum_{\ell=0}^{n}\gamma_{\ell}\lambda^{\ell},$$

(5.2)

and set $p(\lambda)$ and $q(\lambda)$ to be the polynomials defined by (4.8) and (4.9). The goal is now to choose polynomials $p(\lambda)$ and $r(\lambda)$ of degrees $m-1$ and $n$, respectively, such that $r(\lambda)$ is monic and

$$\sup_{\lambda\in[-1,1]}\left|\frac{p(\lambda)}{q(\lambda)}-r(\lambda)\right|=\sup_{\lambda\in[-1,1]}\frac{\left|p(\lambda)-q(\lambda)r(\lambda)\right|}{|q(\lambda)|}$$

(5.3)

is close to zero. We choose $p(\lambda)$ and $r(\lambda)$ such that

$$T_{m+n}(\lambda)/2^{m+n-1}=q(\lambda)r(\lambda)-p(\lambda).$$

(5.4)

This is simply the Euclidean division of the normalized Chebyshev polynomial $T_{m+n}(\lambda)/2^{m+n-1}$ by $q(\lambda)$ with $r(\lambda)$ being identified as the quotient polynomial and $-p(\lambda)$ being identified as the remainder polynomial. Then, assuming (4.2) and using (4.14), we have

$$\sup_{\lambda\in[-1,1]}\left|\frac{p(\lambda)}{q(\lambda)}-r(\lambda)\right|\leq\epsilon_{m}^{(n)},\quad\text{with}\quad\epsilon_{m}^{(n)}=\frac{2}{2^{n}(2d_{\rm min})^{m}}$$

(5.5)

satisfying $\epsilon_{m}^{(n)}\to 0$ as $m\to\infty$, with $n$ being fixed. With constants $\alpha_{k}$ given by (4.10) and constants $\gamma_{\ell}$ being the coefficients of the polynomial $r(\lambda)$, as in (5.2), it follows that

			$\displaystyle\sup_{\mu}\left|\sum_{k=1}^{m}\alpha_{k}f(z_{k})-\sum_{\ell=0}^{n}\gamma_{\ell}M_{\ell}\right|$		(5.6)
			$\displaystyle\quad=\sup_{\mu}\int_{-1}^{1}d\mu(\lambda)\left|\sum_{k=1}^{m}\frac{\alpha_{k}}{\lambda-z_{k}}-\sum_{\ell=0}^{n}\gamma_{\ell}\lambda^{\ell}\right|\leq\epsilon_{m}^{(n)}.$

Supposing any constants $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$ are given, it is easy to get bounds on $v(t)$ given by (3.2) at any time $t$ that incorporate the $n$ known moments $M_{1},M_{2},\ldots,M_{n}$. One introduces an angle $\theta$ and Lagrange multipliers $\gamma_{1},\gamma_{2},\ldots,\gamma_{n}$ and takes the minimum value of

$$\int_{-1}^{1}d\mu(\lambda)\mathop{\rm Re}\nolimits\left[e^{i\theta}\sum_{k=1}^{m}\frac{\alpha_{k}e^{-i\omega_{k}(t-t_{0})}}{\lambda-z(\omega_{k})}+\sum_{\ell=1}^{n}\gamma_{\ell}\lambda^{\ell}\right],$$

(6.1)

as $\mu$ varies over all probability measures supported on $[-1,1]$ with unconstrained moments. The minimum will be achieved by the point masses $\mu=\delta(\lambda-\lambda_{0})$, where $\lambda_{0}$ may take one or more values. Typically we will need to choose the Lagrange multipliers $\gamma_{1},\gamma_{2},\ldots,\gamma_{n}$ (that depend on $\theta$) so that the minimum is achieved at $n$ values $\lambda_{0}=\lambda_{0}^{(\ell)}$, $\ell=1,2,\ldots,n$, and then adjust the measure to be distributed at these points

$$d\mu(\lambda)=\sum_{\ell=1}^{n}\mathrm{w}_{\ell}\delta(\lambda-\lambda_{0}^{(\ell)}),$$

(6.2)

with the non-negative weights $\mathrm{w}_{\ell}$, that sum to $1$, chosen so that the moments take their desired values. Then with this measure we obtain the bound

$$\mathop{\rm Re}\nolimits[e^{i\theta}v(t)]\geq a_{0}\sum_{\ell=1}^{n}\mathrm{w}_{\ell}\mathop{\rm Re}\nolimits\left[e^{i\theta}\sum_{k=1}^{m}\frac{\alpha_{k}e^{-i\omega_{k}(t-t_{0})}}{\lambda_{0}^{(\ell)}-z(\omega_{k})}\right].$$

(6.3)

By varying $\theta$ from $0$ to $2\pi$ we obtain bounds that confine $v(t)$ to a convex region in the complex plane. Of course, if we are only interested in bounding $\mathop{\rm Re}\nolimits[v(t)]$, then it suffices to take $\theta=0$ or $\pi$.

Figure 4 and Figure 5 depict the lower and upper bounds on $\mathrm{Re}[v(t)]$ for two systems ($z(\omega)=2+i/\omega$ in Figure 4, thus mimicking the low frequency dielectric response of a lossy dielectric material, and $z(\omega)=2-2/\omega^{2}$ in Figure 5, thus mimicking the dielectric response of a plasma), when the coefficients $\alpha_{k}$ in (6.3) are chosen such that the bounds are extremely tight at $t_{0}=0$, according to (4.3). For both systems, the bounds on $\mathrm{Re}[v(t)]$ are tighter the higher the amount of pieces of information on the system is incorporated. Notice that the bounds colored in black (the largest ones) correspond to the case where only the zeroth order moment $M_{0}$ of the measure is known but not the value of $a_{0}$: in such a case, as shown by the zoomed graph in the blue box, at $t=0$, the upper bound takes value 1 and the lower bound takes value 0, that are the smallest and the highest values $a_{0}$ can take. On the other hand, when $a_{0}$ is assigned, the value that the corresponding bounds take at $t=0$ is exactly $a_{0}=0.6$, as shown by the zoomed graph in the blue box. The graphs show clearly that, in order to estimate the system parameter $a_{0}$, one has just to measure the response of the system at a specific moment of time $t_{0}$ (if the applied field is carefully chosen).

These are the type of bounds used in [33] to bound the temporal response of two-phase composites in antiplane elasticity. It is not yet clear whether those bounds can be derived from variational principles. In general, in the theory of composites, variational methods have proven to be more powerful than analytic approaches. Variational methods produce tighter bounds that often easily extend to multiphase composites: see the books [13, 49, 37, 2, 48] and references therein. For example, the variational approach gives tighter bounds on the complex permittivity at constant frequency of two-phase lossy composites [24], than the bounds obtained by the analytic approach [35, 8]. It also produces bounds on the complex effective bulk and shear moduli of viscoelastic composites [14, 42]. An exception is bounds that correlate the complex effective dielectric constant at more than two frequencies [36] that have yet to be obtained by a systematic variational approach. Variational bounds in the time domain are available [12, 32], but these are nonlocal in time.

Naturally, if one is interested in the response $v_{0}(t)$ at a given (possibly complex) frequency $\omega_{0}$, the easiest solution is to take an input signal $u_{0}(t)$ at that frequency. However, it might not be easy to experimentally generate a signal at that frequency or it might not be easy to measure the response at that frequency. The problem becomes: find complex constants $\alpha_{1},\alpha_{2},\ldots,\alpha_{m}$ such that

$$\sup_{\lambda\in[-1,1]}\left|\sum_{k=1}^{m}\frac{\alpha_{k}}{\lambda-z_{k}}-\frac{1}{\lambda-z_{0}}\right|\leq\epsilon_{m},$$

(7.1)

with $z_{k}=z(\omega_{k})$, $k=0,1,\dots,m$. Defining the polynomials $p(\lambda)$ and $q(\lambda)$ as in (4.8) and (4.9) one needs to find $p(\lambda)$ of degree $m-1$ such

$$\sup_{\lambda\in[-1,1]}\left|\frac{(\lambda-z_{0})p(\lambda)-q(\lambda)}{(\lambda-z_{0})q(\lambda)}\right|\leq\epsilon_{m}.$$

(7.2)

Proceeding as before we choose

$$(\lambda-z_{0})p(\lambda)=q(\lambda)-b_{m}T_{m-1}(\lambda)\quad\text{with}\,\,b_{m}=q(z_{0})/T_{m-1}(z_{0}),$$

(7.3)

where $b_{m}$ has been chosen so that the polynomial $q(\lambda)-b_{m}T_{m-1}(\lambda)$ has a factor of $(\lambda-z_{0})$. Then the residues of $R(\lambda)=p(\lambda)/q(\lambda)$ are given by

$$\alpha_{k}=-b_{m}\frac{T_{m-1}(z_{k})}{(z_{k}-z_{0})\prod_{j\neq k}(z_{k}-z_{j})}=-\frac{T_{m-1}(z_{k})\prod_{j\neq 0}(z_{0}-z_{j})}{T_{m-1}(z_{0})(z_{k}-z_{0})\prod_{j\neq 0,k}(z_{k}-z_{j})},$$

(7.4)

and

$$\sup_{\lambda\in[-1,1]}|(\lambda-z_{0})p(\lambda)-q(\lambda)|=\sup_{\lambda\in[-1,1]}|b_{m}T_{m-1}(\lambda)|=|b_{m}|,$$

(7.5)

so that (7.1) holds with

$$\epsilon_{m}=\frac{|b_{m}|}{d_{0}\inf_{\lambda\in[-1,1]}|q(\lambda)|},$$

(7.6)

where $d_{0}$ denotes the distance from $z_{0}$ to the interval [-1,1]. Joukowski’s map yields:

$$z_{0}=\frac{1}{2}\left(\zeta_{0}+\frac{1}{\zeta_{0}}\right),\,\,\text{with}\,\,R=|\zeta_{0}|>1,$$

(7.7)

whence

$$T_{m-1}(z_{0})=\frac{1}{2}\left(\zeta_{0}^{m-1}+\frac{1}{\zeta_{0}^{m-1}}\right).$$

(7.8)

Moreover, since $\zeta_{0}$ runs over a circle of radius $R$, we have

$$d_{0}=\inf_{\lambda\in[-1,1]}\frac{|\zeta_{0}^{2}-2\zeta_{0}\lambda+1|}{2}\geq\frac{(R-1)^{2}}{2R},$$

(7.9)

and

$$|T_{m-1}(z_{0})|\geq\frac{1}{2}(R^{m-1}-R^{1-m}),\ \ m\geq 2,$$

(7.10)

implying

$$|b_{m}|\leq\frac{2|q(z_{0})|}{R^{m-1}-R^{1-m}}.$$

(7.11)

All in all, the relevant bound $\epsilon_{m}$ satisfies

$$|\epsilon_{m}|\leq\frac{4R}{(R-1)^{2}}\frac{1}{R^{m-1}-R^{1-m}}\sup_{\lambda\in[-1,1]}\frac{|q(z_{0})|}{|q(\lambda)|}.$$

(7.12)

We obtain an exponential decay $\epsilon_{m}\rightarrow 0$ as $m\rightarrow\infty$ provided the geometry of the loci $z_{1},z_{2},\ldots,z_{m}$ is subject to the following condition: for a positive constant $r<R$, each $z_{j}\in H(r)=H_{1}(r)\cup H_{2}(r)\cup H_{3}(r)$ where

	$\displaystyle H_{1}(r)$	$\displaystyle=$	$\displaystyle\left\{z:\left|\frac{z-z_{0}}{z+1}\right|\leq r,\ \ \mathop{\rm Re}\nolimits z\leq-1\right\},$
	$\displaystyle H_{2}(r)$	$\displaystyle=$	$\displaystyle\left\{z:\left|\frac{z-z_{0}}{\mathop{\rm Im}\nolimits z}\right|\leq r,\ \ \mathop{\rm Re}\nolimits z\in[-1,1]\right\},$
	$\displaystyle H_{3}(r)$	$\displaystyle=$	$\displaystyle\left\{z:\left|\frac{z-z_{0}}{z-1}\right|\leq r,\ \ \mathop{\rm Re}\nolimits z\geq 1\right\}.$		(7.13)

In other words, all of the $z_{j}$ must be close to $z_{0}$ in the precise sense that $z_{j}\in H(r)$. Note that, as shown in Figure 6a, in case $r<1$, $H_{1}$ and $H_{3}$ are sectors of disks, while $H_{2}$ is a portion of an ellipse. For $r\in(1,R)$ these regions are complements of disks/ellipse, containing the point $z_{0}$, as shown in Figure 6c. Some of these regions can be empty, depending on the position of $z_{0}$.