Polynomial methods to construct inputs for uniformly ensemble reachable linear systems

Abstract

This paper is concerned with linear parameter-dependent systems and considers the notion of uniform ensemble reachability. The focus of this work is on constructive methods to compute suitable parameter-independent open-loop inputs for such systems. In contrast to necessary and sufficient conditions for ensemble reachability, computational methods have to distinguish between continuous-time and discrete-time systems. Based on recently derived sufficient conditions and techniques from polynomial approximation, we present two methods for discrete-time single-input linear systems. Moreover, we illustrate that one method can also be applied to certain continuous-time single-input systems.

1 Introduction

Ensemble control is an emerging field in mathematical systems and control theory referring to the task of controlling a large, potentially infinite, number of states, or systems, using a single-input function or a single-feedback controller. Being not a rigorously defined term it subsumes several different scenarios, where in each particular case the problems and the required techniques differ from each other. Ensemble control embraces the scenario of uncertainties in the initial data. This is modeled by a probability distribution of the initial state, and the ensemble control problem leads to a transport problem of density functions and therefore to controllability and observability issues of the Liouville and the Fokker–Plank equation [5, 8, 17, 52, 54]. Another setting tackles uncertainties in the model parameters by addressing controllability problems for parameter-dependent systems, and the goal is to achieve a certain control task by using only a single or a few open-loop inputs which are independent of the (usually unknown) model parameters. The notion ensemble controllability is used for this, cf. [33]. In this context we also mention the notions averaged controllability and moment controllability, cf. [51] and [55], respectively. Motivated by the controllability analysis of the Bloch equation, there is also the concept of asymptotic controllability, cf. e.g., [4, 9].

In this paper we consider families of parameter-dependent systems of the form

$$\begin{aligned} \begin{aligned} \tfrac{ \partial }{ \partial t } x ( t , \theta )&= A ( \theta ) x ( t , \theta ) + B ( \theta ) u ( t ) \\ x _{ t + 1 }( \theta )&= A ( \theta ) x _ t ( \theta ) + B ( \theta ) u _t, \end{aligned} \end{aligned}$$

(1)

where $\theta \in {\textbf{P}}$ is considered as a parameter and the parameter space ${\textbf{P}}\subset {\mathbb {C}}$ is assumed to be a Jordan arc in the complex plane, i.e., ${\textbf{P}}$ is the image of a continuous and bijective function defined on a compact interval. Moreover, we assume that the matrix-valued functions $A: {\textbf{P}}\rightarrow {\mathbb {C}}^{n \times n}$ and $B: {\textbf{P}}\rightarrow {\mathbb {C}}^{n \times m}$ are continuous. Throughout the paper we will use the short notation $(A,B) \in C_{n,n}({\textbf{P}})\times C_{n,m}({\textbf{P}})$. In addition, let $C_n({\textbf{P}})$ denote the space of continuous functions from ${\textbf{P}}$ to ${\mathbb {C}}^n$. As we mainly consider discrete-time systems, we take the initial condition $x(0,\theta )=x_0(\theta )=0$. To express that the solutions to (1) are regarded as functions from the parameter space ${\textbf{P}}$ to ${\mathbb {C}}^n$, we denote them by $\varphi (T,u)(\theta ):=\varphi (T,u,0,\theta )$. Further, we wish to highlight the essential property that the input is independent of the parameters.

The notion of reachability we are considering in this paper is as follows: We say that a pair $(A,B) \in C_{n,n}({\textbf{P}})\times C_{n,m}({\textbf{P}})$ is uniformly ensemble reachable (from zero), if for every $f \in C_n({\textbf{P}})$ and $\varepsilon >0$ there are $T>0$ and $u\in L^1([0,T],{\mathbb {C}}^m)$ or $u = (u_{0},u_{1}, \dots , u_{T-1}), u_{i} \in {\mathbb {C}}^{m}$ such that

$$\begin{aligned} \Vert \varphi (T,u)-f\Vert _{\infty }:= \sup _{\theta \in {\textbf{P}}} \Vert \varphi (T,u,0)(\theta )-f(\theta ) \Vert _{{\mathbb {C}}^n} <\varepsilon . \end{aligned}$$

Note that ensemble reachability is an infinite-dimensional problem (unless ${\textbf{P}}$ is finite) and it is equivalent to approximate reachability of the infinite-dimensional linear system defined by the bounded linear matrix-multiplication operators

$$\begin{aligned} {{{ \mathcal {M}}}}_A:C_n({\textbf{P}}) \rightarrow C_n({\textbf{P}}), \quad {{{\mathcal {M}}}}_A f(\theta )=A(\theta )f(\theta ) \end{aligned}$$

and

$$\begin{aligned} {{{\mathcal {M}}}}_B:{\mathbb {C}}^m \rightarrow C_n({\textbf{P}}),\quad {{{\mathcal {M}}}}_B u =B(\theta )u.\end{aligned}$$

As a consequence of the restriction that the inputs are not allowed to depend on the parameter, it is well known that exact reachability (i.e., $\varepsilon =0$) is never possible provided ${\textbf{P}}$ is not finite, cf. [47, Theorem 3.1.1] and [19, p.244]. We note that in the literature continuous-time systems are mostly considered, and hence the term ensemble controllability is present more frequently. Moreover, it follows from [47, Theorem 3.1.1, Remark 3.1.2] that a continuous-time pair (A, B) is uniformly ensemble reachable if and only if it is completely approximately controllable, i.e., for every $T>0$, for every $\varepsilon >0$ and for every pair $x_0,x_1 \in C_n({\textbf{P}})$ there exists an input $u \in L^1([0,T],{\mathbb {C}}^m)$ such that

$$\begin{aligned} \Vert x_1 - \varphi (T,u,x_0)\Vert _{\infty } < \varepsilon . \end{aligned}$$

We note that this equivalence does not hold for discrete-time systems. Recall that a pair $(A,B) \in C_{n,n}({\textbf{P}})\times C_{n,m}({\textbf{P}})$ is called uniformly ensemble controllable (to zero), if for all $x_0 \in C_n({\textbf{P}})$ and $\varepsilon >0$ there exist $T \ge 0$ and an input u such that

$$\begin{aligned} \Vert \varphi (T,u,x_0)\Vert _{\infty } < \varepsilon . \end{aligned}$$

It is well known that the notions of approximate reachability (from zero) and approximate controllability (to zero) are independent of each other, cf. [18, Lemma 4.1], and none of both implies approximate complete controllability. Also we note that this paper does not cover other notions of controllability, such as asymptotic controllability, cf. [45]. Moreover, recall that the application of [47, Theorem 3.1.1] yields that a pair (A, B) is uniformly ensemble reachable if and only if

$$\begin{aligned} \overline{{\text {span}} \{\theta \mapsto A(\theta )^lb_j(\theta ) \,\, |\, \, j=1,...,m, \, l=0,1,...\} } = C_n({\textbf{P}}), \end{aligned}$$

(2)

where $b_j(\theta )$ denotes the jth column of $B(\theta )$. In operator theory language, uniform ensemble reachability of (A, B) is equivalent to the fact that the matrix multiplication operator ${{{\mathcal {M}}}}_A$ is m-multicyclic and the columns $b_1,...,b_m$ are cyclic vectors. For other characterizations, we refer to [12, Section 6.2] and [16, Section VI.8.a]. As these characterizations have the drawback that they are hard to check in practice, much effort has recently been spent on the derivation of pointwise testable necessary and sufficient conditions, cf. [14, 34].

The ensemble controllability problem for continuous-time parameter-dependent systems is also studied in [1, 2, 6, 32,33,34,35] and [51]. Agrachev and Sarychev consider ensemble controllability for nonlinear drift-free parameter-dependent systems and provide a characterization in terms of Lie brackets. In the same direction, the work of Chen [6] also treats nonlinear systems and considers Lie extensions. We note that these approaches do not apply to the setting in this paper. In [32] a characterization for ensemble controllability for time-varying parameter-dependent linear systems is presented, which is based on the singular value decomposition of the reachability operator. Like condition (2), this condition is hard to check. The recent works [7] and [13] have shown that, if (A, B) satisfy appropriate regularity assumptions, $L^p$- and uniform ensemble controllability, respectively, cannot hold if the parameter space ${\textbf{P}}\subset {\mathbb {R}}^d$, $d \ge 2$ has interior points.

Problem statement   The objective of this paper is as follows. We consider uniformly ensemble reachable pairs $(A,B) \in C_{n,n}({\textbf{P}}) \times C_{n,m}({\textbf{P}})$. Given a desired target family $f \in C_n({\textbf{P}})$ and an $\varepsilon$-neighborhood of it, we treat the problem of how to find a suitable $T>0$ and a suitable input u such that

$$\begin{aligned} \Vert \varphi (T,u) - f\Vert _{\infty } < \varepsilon . \end{aligned}$$

For continuous-time systems this problem has recently been addressed in [31, 43] and [53]. The methods in [43] and [53] are based on techniques from integral equations. Also we mention that the results in [43] require an additional assumption on the target family f so that predefined error bounds can be archived. Besides that, the work [31] uses the adjoint system to set up a continuous-time optimal control problem to derive suitable inputs. However, as the solution formula for discrete-time systems is not given in terms of an integral, the methods just mentioned cannot be applied to discrete-time systems. Moreover, we are not aware of any approaches in the literature to tackle the discrete-time case.

Novelty and main contribution   The scope of the present paper is to fill this gap. More precisely we focus on single-input systems and start with treating the discrete-time case

$$\begin{aligned} x _{ t + 1 }( \theta ) = A ( \theta ) x _ t ( \theta ) + {b} ( \theta ) u _t. \end{aligned}$$

The solution at time $T>0$ for inputs $u=(u_0,...,u_{T-1}) \in {\mathbb {C}}\times \cdots \times {\mathbb {C}}$ is given by

$$\begin{aligned} \varphi (T,u)(\theta ) = \sum _{k=0}^{T-1} \big (A(\theta )\big )^{T-1-k}b(\theta )u_{k}. \end{aligned}$$

Noticing the restriction that the inputs $u_0,...,u_{T-1}$ have to be independent of the parameters, the solution can be expressed in terms of the polynomial

$$\begin{aligned} p: {\mathbb {C}}\rightarrow {\mathbb {C}},\quad p(z)= u_{T-1} + u_{T-2}\,z+\cdots + u_{1}\, z^{T-2}+ u_{0} \,z^{T-1}. \end{aligned}$$

That is, we have

$$\begin{aligned} \varphi (T,u)(\theta ) = p(A(\theta ))b(\theta ). \end{aligned}$$

In this notation, the problem can be stated equivalently as follows. Given $f\in C_n({\textbf{P}})$ and $\varepsilon >0$, how to find a suitable polynomial p such that

$$\begin{aligned} \Vert \varphi (T,u) - f\Vert _{\infty } = \Vert p(A)b -f\Vert _\infty < \varepsilon \, {?} \end{aligned}$$

(3)

This observation is the key to a new strategy to solve the construction problem by using tools from polynomial approximation. In this context, given the target family $f \in C_n({\textbf{P}})$ and the $\varepsilon$-neighborhood of it, the task of finding the amount required inputs and their values is equivalent to derive the degree and the coefficients of an approximating polynomial (in the sense of (3)). Note that this fact is significantly different from the continuous-time case, where the solution formula is given in terms of an integral operator.

The contributions of this paper are made up of two parts. On the one hand, the paper provides a collection of results from (complex) polynomial approximation, for which explicit representations of the corresponding polynomials are available. More precisely, after a short recap of explicit refinements of the classical Weierstrass approximation theorems, the paper provides an explicit construction of the polynomial in Runge’s little theorem, where the degree and the coefficients are determined in terms of the target function and the approximation error. To the best of the author’s knowledge, such an explicit representation was not given in the literature so far. Moreover, the paper gives an explicit treatment of results which are due to Walsh. As these results are required later for the derivation of suitable inputs, the paper also provides a detailed treatment of these results (which is only implicitly contained in the literature). More specifically, we show that the construction of appropriate polynomials can be traced back to Runge’s little theorem by using suitable conformal mappings. We note that the numerical treatment of conformal mappings is a wide area. Hence, this topic is beyond the scope of this paper and we only provide some comments to the relevant literature.

On the other hand, based on the exposition on polynomial approximation, the paper shows various methods to determine appropriate inputs. As there is no explicit characterization of uniform ensemble reachability¹ we take two sufficiency conditions, denoted by (S1) and (S2) in Sect. 3, derived in [14] as a starting point. In this paper we present constructive methods for each condition. Moreover, the paper shows that the methods corresponding to the sufficiency condition (S2) can also be applied to continuous-time single-input systems, while the methods for (S1) are limited to discrete-time systems.

Organization   Let us quickly describe the structure of the paper. Section 2 contains the relevant results from polynomial approximation that are required for the input construction procedures. It is divided into three subsections corresponding to the results of Runge, Walsh and Weierstrass, respectively. Section 3 is devoted to the presentation of the constructive procedures to determine suitable inputs for discrete-time parameter-dependent linear single-input systems. After a short summary of relevant known results, the section has two subsections, where each subsection presents methods corresponding to one of the sufficiency conditions. Section 4 depicts that methods corresponding to sufficiency condition (S2) can also applied to continuous-time systems. Finally, Sect. 5 contains additional perspectives about the presented methods.

Notations   For $x\in {\mathbb {R}}$, let $\lfloor x \rfloor := \max \{m \in {\mathbb {Z}}\,:\, m \le x\}$. For a set $\varOmega$, its interior is denoted by ${\text {int}}\varOmega$ and its closure is denoted by $\overline{\varOmega }$. The distance between two sets $\varOmega ,{\tilde{\varOmega }} \ne \emptyset$ is denoted by

$$\begin{aligned} {\text {d}}\!\,(\varOmega ,{\tilde{\varOmega }}):= \inf \{ \Vert \omega - {\tilde{\omega }}\Vert \,: \, \omega \in \varOmega , \, {\tilde{\omega }} \in {\tilde{\varOmega }}\}. \end{aligned}$$

Moreover, we call $C_n(\varOmega )$ the set of continuous functions $g :\varOmega \rightarrow {\mathbb {C}}^n$ and we say that $g \in C_n^k(\varOmega )$ if g is k-times continuously differentiable. A function $g :\varOmega \rightarrow {\mathbb {C}}^n$ is said to satisfy a Lipschitz condition if there is a $L_g>0$ such that

$$\begin{aligned} \Vert g(z_1) - g(z_2)\Vert _{{\mathbb {C}}^n} \le L_g \, \Vert z_1- z_2\Vert _{\varOmega } \end{aligned}$$

for all $z_1, z_2 \in \varOmega$. The set of functions satisfying a Lipschitz condition is denoted by ${\text {Lip}}_n(\varOmega )$. If $\varOmega \subset {\mathbb {C}}$ is compact and $g :\varOmega \rightarrow {\mathbb {C}}$ is continuous, we define

$$\begin{aligned} M_{g}:= \max _{z \in \varOmega } |g(z)|. \end{aligned}$$

For a function $f:[a,b] \rightarrow {\mathbb {R}}$, the nth Bernstein polynomial is given by

$$\begin{aligned} B_{n,f}(x):= \frac{1}{(b-a)^n}\sum _{k=0}^n {n \atopwithdelims ()k}\, f(a+\tfrac{k}{n}(b-a)) \,(x-a)^k \, (b-x)^{n-k}. \end{aligned}$$

(4)

For a continuous function $f:\partial {\mathbb {D}} \rightarrow {\mathbb {C}}$, its kth Fourier coefficient is denoted by

$$\begin{aligned} {\hat{f}}(k) := \frac{1}{2 \pi } \int _{-\pi }^\pi f(e^{is}) e^{-iks} {\text {d}}\!s. \end{aligned}$$

(5)

Further, we call

$$\begin{aligned} F_{f,n}(z) := \sum _{k=-n+1}^{n-1} \frac{n-|k|}{n} {\hat{f}} (k)\, \, z^k, \end{aligned}$$

(6)

the nth Fejér polynomial.

2 Elements of approximation theory

In this section we discuss results from approximation theory, for which constructive proofs are available. The presented results are due to Bernstein, Runge, Weierstrass and Walsh. Note that the definite result in this context is the famous Mergelyan’s theorem, saying that a continuous function $f :K \rightarrow {\mathbb {C}}$ can be uniformly approximated by polynomials if K is compact, ${\mathbb {C}}\setminus K$ is connected and f is analytic in the interior of K, cf. [24, Chap. III, § 2, Theorem 1]. However, its highly ingenious method of proof is not constructive. Therefore we will present some special cases which can be proved constructively and prepare the ground for the computational methods in Sect. 3. In Sect. 2.1 we recall known explicit versions of the classical Weierstrass approximation theorems. Moreover, in Sect. 2.2 we present an explicit treatment of the approximating polynomial in Runge’s little theorem. In Sect. 2.3 we present a detailed derivation of some of Walsh’s results on complex approximation. Note that these results appeared more than twenty years before Mergelyan’s result.

2.1 Weierstrass Theorems

We start with the Weierstrass approximation theorems and consider the special cases where the continuous function that shall be approximated additionally satisfies a Lipschitz condition. Recall that for a function $f\in C([a,b],{\mathbb {R}})$ the nth Bernstein polynomial is given by

$$\begin{aligned} B_{n,f}(x):= \frac{1}{(b-a)^n}\sum _{k=0}^n {n \atopwithdelims ()k}\, f(a+\tfrac{k}{n}(b-a)) \,(x-a)^k \, (b-x)^{n-k}. \end{aligned}$$

For complex-valued functions $f :[a,b] \rightarrow {\mathbb {C}}$ one can apply the latter to the real and imaginary part. That is, for $f(x) = g(x)+ih(x)$ we consider the complex Bernstein polynomial $B_{n,f}(x) = B_{n,g}(x) + i B_{n,h}(x)$. Next, we recap a result of Gzyl and Palacios that provides an explicit error bound, cf. [27].

Theorem 1

(Gzyl, Palacios (1997)) Let $f:[a,b]\rightarrow {\mathbb {C}}$ satisfy a Lipschitz condition. Then, for $n \ge 3$, the sequence of complex Bernstein polynomials satisfies

$$\begin{aligned} \Vert f - B_{n,f}\Vert _{\infty } \le \sqrt{2} \left( 4\,M_f + \frac{(b-a)\,L_f}{2} \right) \frac{\sqrt{\log n}}{\sqrt{n}}. \end{aligned}$$

Note that, since Gyzl and Palacios consider real-valued functions defined on the unit interval, the constants in the previous statement are adjusted. Similarly, the second or trigonometric Weierstrass approximation theorem can also be proven constructively. It is based on Fejér’s theorem, and we first recall their definition. Let $f:\partial {\mathbb {D}} \rightarrow {\mathbb {C}}$ be a continuous function, then the nth Fejér polynomial is given by

$$\begin{aligned} F_{f,n}(z) := \sum _{k=-n+1}^{n-1} \frac{n-|k|}{n} {\hat{f}} (k)\, \, z^k, \end{aligned}$$

where ${\hat{f}}(k)$ denotes the kth Fourier coefficient of f. From [38, Ch. VIII, Sec. 1, Thm. 1 and Sec. 2, Thm. 3] we recap the following result.

Theorem 2

(Second Weierstrass Theorem) Suppose that $f:\partial {\mathbb {D}} \rightarrow {\mathbb {C}}$ is continuous. Let $(F_n)_{n\in {\mathbb {N}}}$ denote the Fejér polynomials.

• (a) The sequence $(F_n)_{n\in {\mathbb {N}}}$ converges uniformly to f.

• (b) If f satisfies a Lipschitz condition, it holds

  $$\begin{aligned} \sup _{z \in \partial {\mathbb {D}}} |f(z)- F_{f,n}(z)| \le 2 \sqrt{2} \pi \,L_f \cdot \frac{\ln n }{n}, \end{aligned}$$

  where $L_f>0$ denotes the Lipschitz constant.

We note that sharper but less explicit error bounds for the second Weierstrass approximation theorem have been derived in Lorch [36], Nikolski [39] and Telyakovskii [46].

2.2 Runge’s Little Theorem

Another famous result from approximation theory is due to Runge. Before presenting a constructive proof, we have to fix some notation. The presentation of Runge’s little theorem is based on [26] and [42]. Let $\gamma$ denote a closed (piecewise) $C^1$-path. With a slight abuse of notation, we denote its trace also by $\gamma$. Then,

$$\begin{aligned} {\text {ind}}_{\gamma }(z) := \frac{1}{2\pi i} \int _{\gamma } \frac{1}{\xi -z} {\text {d}}\!\xi \,, \quad z\in {\mathbb {C}}\setminus \gamma \end{aligned}$$

denotes the winding number of z with respect to $\gamma$. Moreover, a closed polygon $\tau =[p_1\,p_2\,\cdots \, p_k\, p_1]$ composed of finitely many horizontal or vertical segments $[p_1\,p_2]$,$[p_2\,p_3]$,...,$[p_k\,p_1]$ is called a grid polygon if there exists a not necessarily regular grid $G \subset {\mathbb {C}}$ of horizontal or vertical lines such that all vertices $p_1,...,p_k$ are pairwise distinct adjacent grid points of G. Moreover, let

$$\begin{aligned} {\text {ext}} \gamma := \{ z\in {\mathbb {C}}\setminus \gamma \, | \, {\text {ind}}_{\gamma }(z) = 0\} \end{aligned}$$

and

$$\begin{aligned} {\text {int}} \gamma := \{ z\in {\mathbb {C}}\setminus \gamma \, | \, {\text {ind}}_{\gamma }(z) = 1\} \quad \text { or } \quad {\text {int}} \gamma := \{ z\in {\mathbb {C}}\setminus \gamma \, | \, {\text {ind}}_{\gamma }(z) = -1\} \end{aligned}$$

depending on the orientation of $\gamma$. We will not provide a complete proof. Instead, as we are interested in constructive methods, we will provide an explicit treatment of the degree of the approximating polynomial as well as an explicit representation of its coefficients. To best of the author’s knowledge these are not available in the literature. For some analytical parts that will be left out we refer to [42].

Theorem 3

(Runge’s Little Theorem (1885)) Let $K\subset {\mathbb {C}}$ be compact such that ${\mathbb {C}}\setminus K$ is connected. If there is an open set $\varOmega$ containing K such that f is holomorphic on $\varOmega$, then for every $\varepsilon >0$ there is a polynomial p such that

$$\begin{aligned} \sup _{z\in K} |f(z)-p(z)|< \varepsilon . \end{aligned}$$

Proof

The proof is carried out in four constructive steps. We start with the construction of finitely many grid polygons in $\varOmega {\setminus } K$ so that f can be represented via the Cauchy integral formula for compact sets.

Step 1: Grid polygon construction There are finitely many distinct oriented horizontal or vertical line segments $\tau _1,...,\tau _N$ of length $\delta < \tfrac{1}{\sqrt{2}} {\text {d}}(K,\partial \varOmega )$ so that

$$\begin{aligned} f(z)= \sum _{k=1}^N \frac{1}{2\pi i} \int _{\tau _k} \frac{f(\xi )}{\xi -z} {\text {d}}\!\xi \qquad \forall \, \, z \in K. \end{aligned}$$

Proof of Step 1: Let $\delta \in (0,\tfrac{1}{\sqrt{2}} {\text {d}}(K,\partial \varOmega ))$ and consider a grid consisting of lines that are parallel to the real and imaginary axis so that the distance between two parallel lines is $\delta$. Since K is compact, there are finitely many boxes $Q_1,...,Q_n$ that intersect with K and satisfy

$$\begin{aligned} K \subset \bigcup _{k=1}^n Q_k \subset \varOmega . \end{aligned}$$

Indeed, for $q_l \in Q_l$ it follows from $2\delta < {\text {d}}(K,\partial \varOmega )$ that $B_{\delta }(q_l) \subset \varOmega$. As each box has diameter $\sqrt{2}\delta$ and

$$\begin{aligned} \sqrt{2}\delta< {\text {d}}(K,\partial \varOmega )< {\text {d}}(q_l,\partial \varOmega ){,} \end{aligned}$$

it follows that $Q_l\subset \varOmega$.

The boundaries of the boxes $Q_1,...,Q_n$ consist of line segments. We now pick those boundary segments that are not common for two distinct boxes $Q_l$ and $Q_m$ with $l\ne m$. Let $\tau _1,...,\tau _N$ denote the selected segments. The construction yields that

$$\begin{aligned} \bigcup _{k=1}^N \tau _k \subset \varOmega \setminus K. \end{aligned}$$

Indeed, if line segment $\tau _k$ meets K, there are two boxes of the grid having $\tau _k$ as a common side, a contradiction to the choice of the segments $\tau _1,..., \tau _N$.

Let $z \in K$ be a point that does not lie on any boundary of the boxes $Q_1,...,Q_n$. Then there is exactly one box $Q_l$ containing z and it holds

$$\begin{aligned} f(z)= \frac{1}{2\pi i} \int _{\partial Q_l} \frac{f(\xi )}{\xi -z} {\text {d}}\!\xi \quad \text { and } \int _{\partial Q_k} \frac{f(\xi )}{\xi -z} {\text {d}}\!\xi =0 \quad \forall \, k\ne l . \end{aligned}$$

Since common sides of different boxes of the grid occur twice with different orientation, it follows that

$$\begin{aligned} f(z)= \frac{1}{2\pi i} \sum _{k=1}^n \int _{\partial Q_k} \frac{f(\xi )}{\xi -z} {\text {d}}\!\xi = \frac{1}{2\pi i} \sum _{k=1}^N \int _{\tau _k} \frac{f(\xi )}{\xi z} d\xi . \end{aligned}$$

Moreover, as shown in [42, Chapter 12, § 1.1], the last equality also holds for all $z \in K$ that are on the boundary of some box.   $\Box$

Now, let $\varepsilon >0$, $\delta >0$ and N be fixed. Moreover, let $L>0$ denote the Lipschitz constant of the function $\xi \mapsto \tfrac{f(\xi )}{\xi -z}$ on $\tau$.

Step 2: Rational approximation Let $M:= \lfloor \tfrac{3}{\varepsilon }\tfrac{2\pi }{\delta ^2 \,L\,N}\rfloor$ and divide each line segment in M subsegments of length $\tfrac{\delta }{M}$. Let the subsegments be denoted by $\tau _{k,l}$, for $l=1,...,M$ and $k=1,...,N$. Take distinct points $w^{(k)}_{l} \in \tau _{k,l}$ such that $|w^{(k)}_{l+1} - w^{(k)}_{l}| \le \tfrac{\delta }{M}$ and let

$$\begin{aligned} c^{(k)}_l:= f(w^{(k)}_l) \int _{\tau _{k,l}} {\text {d}}\!\xi \quad \text { for } l=1,...,M \text { and } k=1,...,N. \end{aligned}$$

Then, the rational function

$$\begin{aligned} r(z)= \frac{1}{2\pi i} \sum _{k=1}^{N} \sum _{l=1}^{M} \frac{c^{(k)}_l}{w^{(k)}_l-z} \end{aligned}$$

approximates f uniformly on K, i.e., r satisfies

$$\begin{aligned} \sup _{z\in K} | f(z)-r(z)|< \tfrac{\varepsilon }{3}. \end{aligned}$$

Proof of Step 2: It holds

$$\begin{aligned} \left| \int _{\tau _k} \frac{f(\xi )}{\xi -z}{\text {d}}\!\xi - \sum _{l=1}^{M} \frac{c^{(k)}_l}{w^{(k)}_l-z} \right|&= \left| \sum _{l=1}^M \int _{\tau _{k,l}} \frac{f(\xi )}{\xi -z}{\text {d}}\!\xi - \sum _{l=1}^{M} \int _{\tau _{k,l}} \frac{f(w^{(k)}_l)}{w^{(k)}_l-z} {\text {d}}\!\xi \right| \\&\le \sum _{l=1}^M \left| \int _{\tau _{k,l}} \left( \frac{f(\xi )}{\xi -z} - \frac{f(w^{(k)}_l)}{w^{(k)}_l-z} \right) {\text {d}}\!\xi \right| \\&\le \sum _{l=1}^M \frac{\delta }{M} \max _{\xi \in \tau _{k,l}} \left| \frac{f(\xi )}{\xi -z} - \frac{f(w^{(k)}_l)}{w^{(k)}_l-z} \right| . \end{aligned}$$

As the function $\xi \mapsto \tfrac{f(\xi )}{\xi -z}$ is Lipschitz continuous on $\tau$, one has

$$\begin{aligned} \left| \int _{\tau _k} \frac{f(\xi )}{\xi -z}{\text {d}}\!\xi - \sum _{l=1}^{M} \frac{c^{(k)}_l}{z-w^{(k)}_l} \right|&\le \sum _{l=1}^M \frac{\delta }{M} L\, \max _{\xi \in \tau _{k,l}} \left| \xi - w^{(k)}_l \right| \le \frac{\delta ^2L}{M} \end{aligned}$$

for every $z \in K$. Thus, for all $z\in K$, one has

$$\begin{aligned} |f(z) - r(z)| \le \frac{1}{2\pi } \sum _{k=1}^N \left| \int _{\tau _k} \frac{f(\xi )}{\xi -z} {\text {d}}\!\xi - \sum _{l=1}^{M} \frac{c^{(k)}_l}{w^{(k)}_l-z} \right| \le \frac{N}{2\pi }\cdot \frac{\delta ^2L}{M} \, < \frac{\varepsilon }{3}. \end{aligned}$$

  $\Box$

Step 3: Pole shifting Let $\eta = \max _{z \in K} |z|$ and let $b \in {\mathbb {C}}{\setminus } K$ with $|b|>2\eta$, $\delta _{kl}:= | w^{(k)}_l -b|$ and $\alpha :=\min _{z\in K} |z-b|$. Then, the rational function

$$\begin{aligned} r_b(z)= \frac{1}{2\pi i} \sum _{k=1}^{N} \sum _{l=1}^{M} c^{(k)}_l q_{kl}\left( \frac{1}{z-b} \right) , \end{aligned}$$

defined by the polynomials $q_{kl}(z) =z+ (w^{(k)}_l-b) z^2+ \cdots + (w^{(k)}_l -b)^{m_{kl}} z^{m_{kl}+1}$ of degree

$$\begin{aligned} m_{kl}+1 \ge \frac{\log \left( \frac{\varepsilon }{3} \frac{2\pi (\alpha -\delta _{kl})}{ NM\,|c^{(k)}_l|} \right) }{\log (\delta _{kl}) - \log (\alpha ) } \end{aligned}$$

(7)

approximates r uniformly on K, i.e., $r_b$ satisfies

$$\begin{aligned} \sup _{z\in K} | r(z)-r_b(z)|< \tfrac{\varepsilon }{3}. \end{aligned}$$

Proof of Step 3: Since

$$\begin{aligned} \frac{1}{z-w^{(k)}_l} = \frac{1}{z-b} \cdot \frac{1}{1-\tfrac{w^{(k)}_l-b}{z-b}}= \frac{1}{z-b} \sum _{v=0}^\infty \left( \tfrac{w^{(k)}_l-b}{z-b}\right) ^v \end{aligned}$$

(8)

for all $z \in K$, one has

$$\begin{aligned} |r(z)- r_b(z)|&\le \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} | c^{(k)}_l| \,\left| \frac{1}{w^{(k)}_l-z} - q_{kl}\left( \frac{1}{z-b} \right) \right| \\&= \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} | c^{(k)}_l| \,\left| \frac{1}{w^{(k)}_l-z} - \frac{1}{z-b} \sum _{p=0}^{m_{kl}} \left( \frac{w^{(k)}_l -b}{z-b} \right) ^p \right| \\&\le \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} \, \frac{| c^{(k)}_l| }{|z-b|} \sum _{p=m_{kl}+1}^\infty \left| \frac{w^{(k)}_l -b}{z-b} \right| ^p \\&\le \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} \, \frac{| c^{(k)}_l| }{\alpha } \sum _{p=m_{kl}+1}^\infty \left| \frac{\delta _{kl}}{\alpha } \right| ^p \\&\le \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} \, \frac{| c^{(k)}_l| }{\alpha - \delta _{kl}} \left( \frac{\delta _{kl}}{\alpha }\right) ^{m_{kl}+1}. \end{aligned}$$

The assertion then follows from observing that (7) is equivalent to

$$\begin{aligned} \left( \frac{\delta _{kl}}{\alpha }\right) ^{m_{kl}+1} \le \frac{\varepsilon }{3} \, \frac{2\pi (\alpha -\delta _{kl}) }{ NM\,|c^{(k)}_l|}. \end{aligned}$$

  $\Box$

Step 4: Polynomial approximation The polynomial

$$\begin{aligned} p(z):= \frac{1}{2\pi i} \sum _{k=1}^{N} \sum _{l=1}^{M} c^{(k)}_l \sum _{v=0}^{m_{kl}} (w^{(k)}_l-b)^v p_{kl}^{(v)}(z), \end{aligned}$$

(9)

where

$$\begin{aligned} p_{kl}^{(v)}(z)&:= a^{(klv)}_0 + a^{(klv)}_1z+\cdots + a^{(klv)}_{m^{(v)}_{kl}} z^{m^{(v)}_{kl}+1},\\ a^{(klv)}_s&:= \frac{1}{2\pi i} \int _{\partial B_{r}(0)} \frac{{\text {d}}\!\xi }{\xi ^{s+1} (\xi -b)^{v+1}}, \qquad s=0,1,...,m_{kl}^{(v)}, \quad r:=\tfrac{|b|}{2}, \end{aligned}$$

and

$$\begin{aligned} m^{(v)}_{kl} +1 > \frac{\log \left( 1+\frac{\varepsilon }{3} \cdot \frac{2\pi }{ NM\,r\, |c^{(k)}_l|} \cdot (r-\eta ) \cdot (\delta _{kl}-r)\right) }{\log \delta _{kl} - \log r } ,\qquad \eta := \max _{z \in K}|z|{,} \end{aligned}$$

(10)

approximates $r_b$ uniformly on K, i.e., p satisfies

$$\begin{aligned} \sup _{z\in K} | r_b(z)-p(z)|< \tfrac{\varepsilon }{3}. \end{aligned}$$

Proof of Step 4: By construction, for $z \in K$, we have

$$\begin{aligned} | r_b(z)-p(z)|&\le \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} |c^{(k)}_l| \left| \sum _{v=0}^{m_{kl}} \frac{(w^{(k)}_l-b)^v}{(z-b)^{v+1}} - (w^{(k)}_l-b)^v p_{kl}^{(v)}(z) \right| \\&\le \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} |c^{(k)}_l| \sum _{v=0}^{m_{kl}} | w^{(k)}_l-b|^v \, \left| \frac{1}{(z-b)^{v+1}} - p_{kl}^{(v)}(z) \right| . \end{aligned}$$

Since the function $\frac{1}{(z-b)^{v+1}}$ is holomorphic on $B_r(0)$, its series expansion around zero is given by the coefficients

$$\begin{aligned} a^{(klv)}_\mu := \frac{1}{2\pi i} \int _{\partial B_r(0)} \frac{{\text {d}}\!\xi }{\xi ^{\mu +1} (\xi -b)^{v+1}}, \qquad \mu =0,1,2,3,...., \quad r:=\tfrac{|b|}{2}>1. \end{aligned}$$

That is,

$$\begin{aligned} \frac{1}{(z-b)^{v+1}} = \sum _{\mu =0}^\infty a_\mu ^{(klv)} z^\mu {,} \end{aligned}$$

and one has

$$\begin{aligned} \left| \frac{1}{(z-b)^{v+1}} - p_{kl}^{(v)}(z) \right| \le \sum _{\mu =m_{kl}^{(v)}+1}^\infty |a_\mu ^{(klv)}| |z|^\mu . \end{aligned}$$

By the Cauchy estimates [26, Theorem 3.4.1], one has

$$\begin{aligned} |a^{(klv)}_\mu | \le \frac{1}{r^{\mu }} \max _{|\xi |=r} \frac{1}{|\xi -b|^{v+1}} = \frac{1}{r^{\mu }} \cdot \frac{1}{\min _{|\xi |=r} |\xi -b|^{v+1}} = \frac{1}{r^{\mu }} \cdot \frac{1}{ r^{v+1}} . \end{aligned}$$

Consequently, with $| w^{(k)}_l-b|^v = (\delta _{kl})^v$ for $z \in K$, we have

$$\begin{aligned} | r_b(z)-p(z)|&\le \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} |c^{(k)}_l| \, \sum _{v=0}^{m_{kl}} \, \frac{(\delta _{kl})^v}{r^{v+1}} \sum _{\mu =m_{kl}^{(v)}+1}^\infty \left( \tfrac{|z|}{r}\right) ^\mu \\&\le \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} |c^{(k)}_l| \, \sum _{v=0}^{m_{kl}} \, \frac{(\delta _{kl})^v}{r^{v+1}} \sum _{\mu =m_{kl}^{(v)}+1}^\infty \left( \tfrac{\eta }{r}\right) ^\mu \\&= \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} |c^{(k)}_l|\, \frac{\left( \tfrac{\eta }{r}\right) ^{m_{kl}^{(v)}+1}}{r-\eta } \sum _{v=0}^{m_{kl}} \, \left( \frac{\delta _{kl}}{r} \right) ^{v} \\&= \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} |c^{(k)}_l|\, \frac{\left( \tfrac{\eta }{r}\right) ^{m_{kl}^{(v)}+1}}{r-\eta } \, \left( \frac{ 1- \left( \tfrac{\delta _{kl}}{r}\right) ^{m_{kl}+1} }{1-\frac{\delta _{kl}}{r} } \right) \\&< \frac{1}{2\pi } \sum _{k=1}^{N} \sum _{l=1}^{M} |c^{(k)}_l|\, \frac{r}{(r-\eta )(r-\delta _{kl})} \, \left( 1- \left( \tfrac{\delta _{kl}}{r}\right) ^{m_{kl}+1} \right) , \end{aligned}$$

where in the last inequality we used that $\left( \tfrac{\eta }{r}\right) ^{m_{kl}^{(v)}+1}<1$ since $\eta < r$. By construction, it holds $\delta _{kl}>r$ and thus (10) is equivalent to

$$\begin{aligned} 1 - \left( \frac{\delta _{kl}}{r}\right) ^{m_{kl}^{(v)}+1} < \frac{\varepsilon }{3} \cdot \frac{2\pi }{NM\,r\,|c^{(k)}_l|} \cdot (r-\eta ) \cdot (r-\delta _{kl}) . \end{aligned}$$

This completes step 4.   $\Box$

The above construction steps yield for every $z\in K$

$$\begin{aligned} |f(z) - p(z)| \le |f(z) - r(z)| + |r(z) - r_b(z)| + |r_b(z) - p(z)| < \varepsilon , \end{aligned}$$

which completes the proof. $\blacksquare$

Remark 1

• (a) The set $\bigcup _{k=1}^N \tau _k$ is determined only by the sets $\varOmega$ and K, and the grid constructed in Step 1. The set $\bigcup _{k=1}^N \tau _k$ also contains all the poles of r and is independent of the quality of approximation.

• (b) As $\varepsilon >0$ is decreasing, the number of poles of the approximating rational function is increasing. The location of the poles is not affected, i.e., the poles will still be located on $\bigcup _{k=1}^N \tau _k$.

• (c) It is shown in [42, Chapter 12, § 4] that the line segments $\tau _1,...,\tau _N$ define a cycle $\tau$ consisting of finitely many grid polygons. In general, it is not possible to find a single closed path $\varGamma$ such that

  $$\begin{aligned} f(z)= \frac{1}{2\pi i} \int _{\varGamma } \frac{f(\xi )}{\xi -z} {\text {d}}\!\xi \qquad \forall \, \, z \in K. \end{aligned}$$

  A counterexample is given [42, Chapter 12, § 1]. However, under suitable assumptions on the compact set K, the line segments $\tau _1,...,\tau _N$ yield a single closed grid polygon, cf. [14, Theorem 9 (b)].

2.3 Walsh’s Theorems

In this subsection we shall prove Walsh’s refinement of Runge’s theorem. We note that Walsh’s result is still considerably short of Mergelyan’s definite result on complex approximation, cf. [24, Chapter III §2 A]. We note that the proof of Mergelyan’s result is not constructive, which is true for almost all results in polynomial approximation. Recall that, a Jordan arc is a homeomorphic image of a compact interval. The result is as follows, cf. [50, Bemerkung $6^{\circ }$].

Theorem 4

(Walsh [50]) Let $\gamma$ be a Jordan arc. If f is continuous on $\gamma$, it can be uniformly approximated on $\gamma$ by a polynomial in z.

A proof of Theorem 4 is given at the end of this section. We emphasize that the following exposition is more explicit in terms of constructing the polynomial than the ones in the literature. The significance of the subsequent exposition is that the construction of the polynomial p is traced back to Runge’s result by using suitable conformal mappings. Numerical procedures for conformal mappings are an active research area and beyond the scope of this paper. At the end of this section we provide more detailed comments and references.

Before we come to the proof, we need some preparation. To do so, we have to consider sequences of Jordan curves, or more precisely we provide a notion of convergence for sequences of Jordan curves. Recall that, a Jordan curve is a homeomorphic image (within ${\mathbb {C}}$) of the unit circle $\partial {\mathbb {D}}$.

Following Courant [11], a sequence $(\varGamma _n)_{n \in {\mathbb {N}}} \subset {\mathbb {C}}$ of Jordan curves is said to converge to a Jordan curve $\varGamma$ if (i) any accumulation point of any sequence $({z}_n)_{n \in {\mathbb {N}}}$ with ${z}_n \in {\text {int}}\varGamma _n$ lies on $\varGamma$ and (ii) for every $\varepsilon >0$ there are $n(\varepsilon ) \in {\mathbb {N}}$ and $\delta (\varepsilon )>0$ with $\lim _{\varepsilon \rightarrow 0} \delta (\varepsilon ) = 0$ such that for every ${z} \in \varGamma$ and for every ${z}_1,{z}_2 \in {\text {int}}\varGamma _n \cap B_{\varepsilon }(p)$, $n > n(\varepsilon )$ there is a polygon $\tau _n$ in $\varGamma _n$ connecting ${z}_1$ and ${z}_2$ so that ${\text {diam}} \tau _n:= \sup \{|{v}-w|\, {:\, v}, w \in \tau _n\}< \delta (\varepsilon )$.

Let $\varOmega$ be a proper subdomain of ${\mathbb {C}}$. Recall that a function $f:\varOmega \rightarrow {\mathbb {D}}$ is called conformal, if it is holomorphic, one-to-one and onto. The following result is due to Courant, cf. [10, 11].

Theorem 5

(Courant (1914)) Let $\varOmega$ and $(\varOmega _n)_{n\in {\mathbb {N}}}$ be domains containing the origin such that the boundaries $\varGamma =\partial \varOmega$ and $\varGamma _n=\partial \varOmega _n$ are Jordan curves. Let $\psi :\overline{{\mathbb {D}}} \rightarrow \overline{\varOmega }$ and $\psi _n:\overline{{\mathbb {D}}} \rightarrow \overline{\varOmega _n}$, $n\in {\mathbb {N}}$, denote continuous functions that are one-to-one, onto, and conformal on ${\mathbb {D}}$ such that $\psi (0)=\psi _n(0)=0$ and $\psi '(0)>0$ and $\psi '_n(0)>0$. Then $(\psi _n)_{n \in {\mathbb {N}}}$ converges uniformly on $\overline{{\mathbb {D}}}$ to $\psi$ if and only if $(\varGamma _n)_{n\in {\mathbb {N}}}$ converges to $\varGamma$.

Using a different method, Courant’s theorem was also derived by Radó [41]. We note that Radó used the Fréchet-distance to characterize convergence for sequences of Jordan curves. That these notions are in fact equivalent was shown by Markushevich. A nice and in-depth exposition for this, as well as extensions to more general domains can be found in [21].

For a given Jordan curve $\varGamma$, a sequence of Jordan curves $(\varGamma _n)_{n \in {\mathbb {N}}}$ that converges to $\varGamma$ from outside, i.e., $\varGamma _n \subset {\mathbb {C}}\setminus \overline{ {\text {int}}\varGamma }$, can be constructed via the grid construction in Step 1 in the proof Runge’s theorem. Indeed, as illustrated in Step 1 in the proof Runge’s theorem, it is sufficient to choose a grid consisting of horizontal and vertical lines defined by

$$\begin{aligned} x = \frac{k}{2^n} \qquad y= \frac{l}{2^n}, \quad k,l \in {\mathbb {Z}} \end{aligned}$$

and define $\varGamma _n$ as the boundary of the boxes given by the vertices

$$\begin{aligned} \left( \frac{k}{2^n}, \frac{l}{2^n} \right) , \, \, \left( \frac{k+1}{2^n}, \frac{l}{2^n} \right) , \, \, \left( \frac{k+1}{2^n}, \frac{l+1}{2^n} \right) , \, \,\left( \frac{k}{2^n}, \frac{l+1}{2^n} \right) \end{aligned}$$

that intersect with $\overline{ {\text {int}}\varGamma }$. Next we recap a proof for a refinement of Runge’s theorem which is also due to Walsh, cf. [49, Satz].

Theorem 6

(Walsh [49]) Let $\varGamma$ be a Jordan curve and suppose that f is holomorphic on ${\text {int}}\varGamma$ and continuous on $\overline{ {\text {int}}\varGamma }$. Then, f can be uniformly approximated on $\overline{ {\text {int}}\varGamma }$ by polynomials. That is, for every $\varepsilon >0$ there is a polynomial p such that

$$\begin{aligned} \sup _{z\in \overline{ {\text {int}}\varGamma } } | f (z)- p(z)| < \varepsilon . \end{aligned}$$

Proof

Let $\varepsilon >0$ and let w.l.o.g. $z=0 \in {\text {int}}\varGamma$. Using the grid construction as in the proof of Runge’s theorem there is a sequence $(\varGamma _n)_{n \in {\mathbb {N}}}$ of grid polygons such that $\varGamma _n \subset {\mathbb {C}}{\setminus } \overline{ {\text {int}}\varGamma }$ and $(\varGamma _n)_{n \in {\mathbb {N}}}$ converges to $\varGamma$. Then we consider the conformal mappings

$$\begin{aligned} \varPhi&:{ {\text {int}}\varGamma } \rightarrow {\mathbb {D}}, \quad \,\,\,\, \, \varPhi (0)=0, \, \varPhi '(0)>0 \\ \varPhi _n&:{ {\text {int}}\varGamma _n} \rightarrow {\mathbb {D}}, \quad \varPhi _n(0)=0, \, \varPhi _n'(0)>0 . \end{aligned}$$

Step 1: (Approximation of f by a sequence of holomorphic functions on $\overline{ {\text {int}}\varGamma }$).

We consider the holomorphic functions

$$\begin{aligned} g_n&{:=} f \circ \varPhi ^{-1}\circ \varPhi _n :{ {\text {int}}\varGamma _n} \rightarrow {\mathbb {C}}\end{aligned}$$

and shall show that for every $\varepsilon >0$ there is an $N\in {\mathbb {N}}$ such that

$$\begin{aligned} |g_N(z) - f(z)| < \tfrac{\varepsilon }{2} \qquad \text { for all } \quad z \in \overline{ {\text {int}}\varGamma }. \end{aligned}$$

(11)

Indeed, for $n \in {\mathbb {N}}$ we consider the Jordan curve $\gamma _n:=\varPhi _n(\varGamma ) \subset {\mathbb {D}}$ and note that for every $z \in \overline{ {\text {int}}\varGamma }$ there exists an $w\in \overline{ {\text {int}}\gamma _n}$ so that $z=\varPhi _n^{-1}(w)$. Thus, together with the maximum modulus theorem [26, Corollary 5.3.4], it is sufficient to show that there is a $K=K(\varepsilon ) \in {\mathbb {N}}$ such that

$$\begin{aligned} \sup _{w \in \gamma _K } |f( \varPhi ^{-1}(w))- f( \varPhi _K^{-1}(w))| < \tfrac{\varepsilon }{2}. \end{aligned}$$

Note that, $\overline{ {\text {int}}\varGamma }$ is compact and f is continuous, and so f is uniformly continuous on $\overline{ {\text {int}}\varGamma }$. Thus there is a $\delta (\varepsilon )>0$ such that $|f(z)-f(w)|< \tfrac{\varepsilon }{2}$ for all $z,w \in \overline{ {\text {int}}\varGamma }$ with $|z-w|<\delta (\varepsilon )$. Since the Jordan curves $(\varGamma _n)_{n\in {\mathbb {N}}}$ converge to $\varGamma$, we can apply Courant’s Theorem 5 to the sequence $(\varPhi ^{-1}_n)_{n\in {\mathbb {N}}}$ and conclude that there is an $N=N(\varepsilon ) \in {\mathbb {N}}$ such that

$$\begin{aligned} |\varPhi ^{-1}(w) - \varPhi ^{-1}_N(w)| < \delta (\varepsilon ) \qquad \text { for all } \quad w \in \gamma _N \subset {\mathbb {D}}, \end{aligned}$$

which shows (11).

Step 2: (Approximation of$g_N$ via Runge’s theorem). The function

$$\begin{aligned} g_N&= f \circ \varPhi ^{-1}\circ \varPhi _N :\overline{ {\text {int}}\varGamma } \rightarrow {\mathbb {C}}\end{aligned}$$

is holomorphic on the compact set $\overline{ {\text {int}}\varGamma }$. Then, by Runge’s Theorem 3, for $\tfrac{\varepsilon }{2}>0$ there is a polynomial p such that

$$\begin{aligned} \sup _{z \in \overline{ {\text {int}}\varGamma } } |g_N(z)- p(z)| < \tfrac{\varepsilon }{2}. \end{aligned}$$

Consequently, for every $z \in \overline{ {\text {int}}\varGamma }$ it follows that

$$\begin{aligned} |f(z)-p(z)|\le |f(z)-g_N(z)| + |g_N(z)- p(z)| < \varepsilon . \end{aligned}$$

$\blacksquare$

We will use the latter statement to prove Walsh’s extension of the second Weierstrass theorem to arbitrary Jordan curves, cf. [50, Satz I].

Theorem 7

(Walsh [50]) Let $\varGamma$ be a Jordan curve containing the origin in its interior and let the function f be continuous on $\varGamma$. Then f can be uniformly approximated on $\varGamma$ by polynomials in z and $\tfrac{1}{z}$.

Proof

Let $\varepsilon >0$. By Caratheodory’s theorem [26, Theorem 13.2.3], there is a continuous and bijective map $\varPhi :\overline{ {\text {int}}\varGamma } \rightarrow \overline{{\mathbb {D}}}$ so that $\varPhi (0)=0$ which is conformal in ${\mathbb {D}}$. Then, the mapping

$$\begin{aligned} g:=f\circ \varPhi ^{-1} :\partial {\mathbb {D}}\rightarrow {\mathbb {C}}\end{aligned}$$

is continuous. Consider the sequence of Fejér polynomials $(F_{g,n})_{n\in {\mathbb {N}}}:\partial {\mathbb {D}}\rightarrow {\mathbb {C}}$ defined by

$$\begin{aligned} F_{g,n}(w) := \sum _{k=1}^{n-1} \frac{n-|k|}{n} {\hat{g}} (-k)\, \, (\tfrac{1}{w})^{k} + \sum _{k=0}^{n-1} \frac{n-|k|}{n} {\hat{g}} (k)\, \, w^{k}, \end{aligned}$$

where the ${\hat{g}}(k)$ denote the Fourier coefficients of g, cf. (5). Since for every $z \in \varGamma$ there is a unique $w \in \partial {\mathbb {D}}$ such that $z= \varPhi ^{-1}(w)$, by Theorem 2 (a), there is a $N=N(\varepsilon ) \in {\mathbb {N}}$ such that

$$\begin{aligned} \sup _{z \in \varGamma }|f(z)- F_{g,N}(\varPhi (z))| = \sup _{w\in \partial {\mathbb {D}}} |g(w) -F_{g,N}(w)| < \frac{\varepsilon }{3}. \end{aligned}$$

To show the claim, we shall prove that on $\varGamma$, the function $\varPhi$ can be uniformly approximated on $\varGamma$ by polynomials in z and the function $\tfrac{1}{\varPhi }$ can be uniformly approximated by polynomials in z and $\tfrac{1}{z}$.

Let $\delta >0$; we will specify it later, cf. (13). By Theorem 6, there is a polynomial $p_1$ such that

$$\begin{aligned} | {\varPhi (z)} - p_1(z)| < \delta \qquad \forall \, \, z \in \overline{ {\text {int}}\varGamma }. \end{aligned}$$

Next, we show that there is a polynomial $p_2$ such that

$$\begin{aligned} \left| \frac{1}{\varPhi (z)}- \frac{a_{-1}}{z} + p_2(z)\right| < \delta \qquad \text { for all } \, \, z \in \varGamma , \end{aligned}$$

where $a_{-1}$ is the residue of $\tfrac{1}{\varPhi }$ at zero, i.e.,

$$\begin{aligned} a_{-1} = \frac{1}{2\pi i} \int _{\partial B_r(0)} \frac{1}{\varPhi (\xi )} {\text {d}}\!\xi , \qquad r >0 \, \, \text {such that } \overline{B_r(0)} \subset {\text {int}}\varGamma . \end{aligned}$$

To see this, recall that $\varPhi$ is conformal on ${\text {int}}\varGamma$ with $\varPhi (0)=0$. Hence, the function $\tfrac{1}{\varPhi }$ is meromorphic on ${\text {int}}\varGamma$ with a simple pole at 0. Moreover, the function $z \mapsto \tfrac{1}{\varPhi (z)} - \tfrac{a_{-1}}{z}$ has a removable singularity at zero. Let $\Psi : {\text {int}}\varGamma \rightarrow {\mathbb {C}}$ denote its holomorphic extension, i.e.,

$$\begin{aligned} \Psi (z) = {\left\{ \begin{array}{ll} \frac{1}{\varPhi (z)} - \frac{a_{-1}}{z} &{} z \in {\text {int}}\varGamma \setminus \{0\}\\ a_0= \frac{1}{2\pi i} \int _{\partial B_r(0)} \frac{1}{\xi \varPhi (\xi )} {\text {d}}\!\xi &{} z=0, \end{array}\right. } \end{aligned}$$

(12)

which is continuous on $\overline{ {\text {int}}\varGamma }$. Then, by Theorem 6, there is a polynomial $p_2$ such that

$$\begin{aligned} \left| \Psi (z)- p_2(z)\right| < \delta \qquad \text { for all } \, \, z \in \overline{ {\text {int}}\varGamma }. \end{aligned}$$

Next, we shall show that, if we choose $\delta >0$ such that

$$\begin{aligned} \sum _{k=0}^{N-1} \left( (1+\delta )^k -1 \right) < \frac{\varepsilon }{3\,c}, \qquad c := \max _{k=-N+1,...,N-1} {\Big |}\tfrac{N-|k|}{N} {\hat{g}}(k){\Big |} \end{aligned}$$

(13)

one has

$$\begin{aligned} \left| \sum _{k=0}^{N-1} \frac{N-|k|}{N} {\hat{g}} (k)\, \, \varPhi (z)^k - \sum _{k=0}^{N-1} \frac{N-|k|}{N} {\hat{g}} (k)\, \, p_1(z)^k\right| < \frac{\varepsilon }{3} \qquad \forall \, \, z\in \varGamma . \end{aligned}$$

(14)

To see this, note that $|\varPhi (z)| =1$ for all $z \in \varGamma$ and, thus,

$$\begin{aligned} |p_1(z)| < 1+ \delta \qquad \text { for all } z \in \varGamma . \end{aligned}$$

Moreover, we recall that for all $v,w \in {\mathbb {C}}$ it holds

$$\begin{aligned} |v^n-w^n| \le |v-w| \sum _{k=0}^{n-1} |v|^k \, |w|^{n-1-k}. \end{aligned}$$

The combination of the latter yields that for all $z\in \varGamma$ one has

$$\begin{aligned}&\left| \sum _{k=0}^{N-1} \frac{N-|k|}{N} {\hat{g}} (k)\, \, \varPhi (z)^k - \sum _{k=0}^{N-1} \frac{N-|k|}{N} {\hat{g}} (k)\, \, p_1(z)^k \right| \\&\quad \le \sum _{k=0}^{N-1} \left| \frac{N-|k|}{N} {\hat{g}} (k)\right| \left| \varPhi (z)^k - p_1(z)^k\right| \\&\quad \le \sum _{k=0}^{N-1} \left| \frac{N-|k|}{N} {\hat{g}} (k)\right| \left| \varPhi (z) - p_1(z)\right| \,\, \sum _{l=0}^{k-1} \,| p_1(z)|^l \le c \sum _{k=0}^{N-1} \left\{ (1+\delta )^k -1 \right\} < \frac{\varepsilon }{3}. \end{aligned}$$

Similarly, for all $z\in \varGamma$ we get

$$\begin{aligned} \left| \sum _{k=1}^{N-1} \frac{N-|k|}{N} {\hat{g}} (-k)\, \, (\tfrac{1}{\varPhi (z)})^k - \sum _{k=1}^{N-1} \frac{N-|k|}{N} {\hat{g}} (-k)\, \, \left( \tfrac{a_{-1}}{z} + p_2(z)\right) ^k \right| < \frac{\varepsilon }{3}. \end{aligned}$$

Finally, defining the polynomial

$$\begin{aligned} P(z,\tfrac{1}{z}) = \sum _{k=1}^{N-1} \frac{N-|k|}{N} {\hat{g}} (-k)\, \, \left( \tfrac{a_{-1}}{z} + p_2(z)\right) ^k + \sum _{k=0}^{N-1} \frac{N-|k|}{N} {\hat{g}} (k)\, p_1(z)^k{,} \end{aligned}$$

we derive the following estimate

$$\begin{aligned} \left| f(z) - P(z,\tfrac{1}{z}) \right|&\le \left| f(z) - F_{g,N}(\varPhi (z)) \right| \\&\quad + \left| \sum _{k=1}^{N-1} \frac{N-|k|}{N} {\hat{g}} (-k)\,\left( \, \left( \tfrac{1}{\varPhi (z)} \right) ^k - \, \, \left( \tfrac{a_{-1}}{z} + p_2(z)\right) ^k \right) \right| \\&\quad + \left| \sum _{k=0}^{N-1} \frac{N-|k|}{N} {\hat{g}} (k)\left( \varPhi (z)^k - p_1(z)^k \right) \right| < 3\, \tfrac{\varepsilon }{3} = \varepsilon . \end{aligned}$$

This shows the assertion. $\blacksquare$

Proof of Theorem 4 Without loss of generality we assume that $\gamma$ does not contain the origin². Then, the Jordan arc can be extended to a Jordan curve $\varGamma$ so that $0 \in {\text {int}}\varGamma$. Also the function f can be continuously extended to $\varGamma$. This can be achieved by defining

$$\begin{aligned} f(z)= w_1 + (w_2-w_1)\frac{z-z_1}{z_2-z_1} \qquad \forall \, \, z \in \varGamma \setminus \gamma , \end{aligned}$$

(15)

where $w_1$ and $w_2$ are the values of f at the end-points $z_1$ and $z_2$ of $\gamma$, respectively. By Caratheodory’s theorem [26, Theorem 13.2.3] there is a continuous and bijective map $\varPhi :\overline{ {\text {int}}\varGamma } \rightarrow \overline{{\mathbb {D}}}$ so that $\varPhi (0)=0$ which is conformal on ${\text {int}}\varGamma$.

Let $\varepsilon >0$. By Theorem 7, for the continuous function $f :\varGamma \rightarrow {\mathbb {C}}$ there are $N=N(\varepsilon ) \in {\mathbb {N}}$ and $a_{-1} \in {\mathbb {C}}$ as well as polynomials $p_1,p_2$, such that the polynomial

$$\begin{aligned} P_{f\circ \varPhi ^{-1},N}(z,\tfrac{1}{z})&:= \sum _{k=1}^{N-1} c_{-k}\, \, \left( \tfrac{a_{-1}}{z} + p_2(z)\right) ^k + \sum _{k=0}^{N-1} c_k\, p_1(z)^k \end{aligned}$$

in z and $\tfrac{1}{z}$ satisfies

$$\begin{aligned} \sup _{z \in {\varGamma }} |f(z) -P_{f\circ \varPhi ^{-1},N}(z,\tfrac{1}{z}) |< \frac{\varepsilon }{2}. \end{aligned}$$

Here, the coefficients $c_{-N+1},...,c_{N-1}$ are given by the Fourier coefficients corresponding to the function $f \circ \varPhi ^{-1}:\partial {\mathbb {D}}\rightarrow {\mathbb {C}}$, cf. (5), i.e.,

$$\begin{aligned} c_k:= \tfrac{N-|k|}{N} \widehat{f\circ \varPhi ^{-1}} (k) , \qquad k \in \{-N+1,...,N-1\}. \end{aligned}$$

Furthermore, on $\gamma$ the function $z \mapsto \tfrac{1}{z}$ can be uniformly approximated by polynomials in z. To see this, let $\gamma _{0,\infty }$ be a Jordan arc connecting the origin and $\infty$ such that $\gamma _{0,\infty } \cap \gamma = \emptyset$. Then, $\varOmega := {\mathbb {C}}\setminus \gamma _{0,\infty }$ is simply connected, $\infty \not \in {\text {int}}\varOmega$ and the function $z \mapsto \tfrac{1}{z}$ is holomorphic on ${\text {int}}\varOmega$. Let

$$\begin{aligned} \nu := \max _{z \in \gamma } \left| \tfrac{1}{z} \right| , \, \, \mu := \max _{z \in \gamma } \left| p_2(z) \right| \quad \text {and} \quad c:= \max _{k=-N+1,....,-1} \left| c_k \right| . \end{aligned}$$

Moreover, choose $\delta >0$ such that

$$\begin{aligned} \delta \cdot \frac{1-(2\delta + |a_{-1}|\nu +\mu )^{N-1}}{1-(2\delta + |a_{-1}|\nu +\mu )} < \frac{\varepsilon }{2 c(N-1)}. \end{aligned}$$

(16)

By Runge’s Theorem 3, there is a polynomial $p_3$ in z such that

$$\begin{aligned} \sup _{z \in \gamma } |\tfrac{a_{-1}}{z} -p_3(z)|< \delta . \end{aligned}$$

Next, we define the polynomial

$$\begin{aligned} p(z) := \sum _{k=0}^{N-1} c_k\,p_1(z)^k + \sum _{k=1}^{N-1} c_{-k} \, \big ( p_3(z) + p_2(z)\big )^k. \end{aligned}$$

(17)

Similar to the proof of Theorem 7, we get

$$\begin{aligned}&\left| \sum _{k=1}^{N-1} c_{-k}\,\big ( \tfrac{a_{-1}}{z}+ p_2(z)\big )^k - \sum _{k=1}^{N-1} c_{-k}\, \big (p_3(z) + p_2(z)\big )^k \right| \\ {}&\quad \le c \sum _{k=1}^{N-1} \left| ( \tfrac{a_{-1}}{z}+ p_2(z))^k - (p_3(z) + p_2(z))^k \right| \\ {}&\quad \le c\, \delta \sum _{k=1}^{N-1} \sum _{l=0}^{k-1} \left| \tfrac{a_{-1}}{z}+ p_2(z)\right| ^l \left| p_3(z) + p_2(z) \right| ^{k-1-l}. \end{aligned}$$

Using the triangle inequality, for every $z \in \gamma$, one has

$$\begin{aligned} | p_3(z)|< \delta + |\tfrac{a_{-1}}{z}| \le \delta + |a_{-1}|\nu \end{aligned}$$

and $| p_3(z) + p_2(z) |< \delta + |a_{-1}|\nu + \mu$ as well as

$$\begin{aligned} |\tfrac{a_{-1}}{z} + p_2(z) | \le |\tfrac{a_{-1}}{z} - p_3(z)| + | p_3(z) - p_2(z) | < 2 \delta + |a_{-1}|\nu + \mu . \end{aligned}$$

Thus, by (16) we have

$$\begin{aligned} \left| \sum _{k=1}^{N-1} c_{-k}\,\big ( \tfrac{a_{-1}}{z}+ p_2(z)\big )^k - \sum _{k=1}^{N-1} c_{-k}\, \big (p_3(z) + p_2(z)\big )^k \right| \\ \le c\, \delta \, (N-1) \sum _{k=0}^{N-2} \left( 2 \delta + |a_{-1}|\nu + \mu \right) ^{k} < \frac{\varepsilon }{2}. \end{aligned}$$

Consequently, it holds

$$\begin{aligned}&\sup _{z \in \gamma } |f(z) -p(z)| \\ {}&\quad \le \sup _{z \in \varGamma } |f(z) -p(z)| \le \sup _{z \in \varGamma } |f(z) -P_{f\circ \varPhi ^{-1},N}(z,\tfrac{1}{z})| \\ {}&\qquad + \sup _{z \in \varGamma } \left| \sum _{k=1}^{N-1} c_{-k}\,\big ( \tfrac{a_{-1}}{z}+ p_2(z)\big )^k - \sum _{k=1}^{N-1} c_{-k}\, \big (p_3(z) + p_2(z)\big )^k \right| <\varepsilon . \end{aligned}$$

This shows the assertion. $\blacksquare$

We close this section with some comments on the construction of conformal mappings and their numerical computation. There are different approaches available in the literature. Peter Henrici presents a method to determine conformal mappings by solving a Dirichlet boundary problem. The solution of it together with its harmonic conjugate determines the conformal map, cf. [29, Theorem 16.5a]. Moreover, Henrici describes numerical techniques to solve the Dirichlet problem by means of numerical methods to solve Symm’s integral equation of first kind, cf. [29, § 16.6 V]. Besides that, Tobin. A. Driscoll and Llyod N. Trefethen depict procedures, as suggested by Peter Henrici, to compute conformal maps by using the Schwarz-Christoffel mapping, cf. [15]. Another, recent approach to compute conformal maps is the Zipper algorithm [37]. For earlier approaches to this area we refer to [22, 23]. Furthermore, approximation techniques for conformal mappings based on Bergman kernels are outlined in [24] and [26]. Finally, we mention the monograph [40] which contains a profound literature review and many more references.

3 Computational methods for discrete-time single-input systems

In this section we consider single-input pairs $(A,b)\in C_{n,n}({\textbf{P}}) \times C_n({\textbf{P}})$, where the parameter space ${\textbf{P}}$ is a Jordan arc. In particular, for given $f\in C_n({\textbf{P}})$ and $\varepsilon >0$, we tackle the central problem of how to determine suitable inputs such that the origin is steered into the $\varepsilon$-neighborhood of f. The basic tools for the constructions will be the results of Bernstein, Runge, Weierstrass and Walsh from the previous section. Indeed, as pointed out in introduction, the solution to

$$\begin{aligned} \begin{aligned} x _{ t + 1 }( \theta ) = A ( \theta ) x _ t ( \theta ) + b ( \theta ) u _t \end{aligned} \end{aligned}$$

can be represented by a polynomial $p(z)= u_0 z^t + u_1z^{t-1} + \cdots +u_{t-1}$, i.e.,

$$\begin{aligned} \varphi (t,u,0) = p(A)b. \end{aligned}$$

Hence, given a target function $f \in C_n({\textbf{P}})$ and $\varepsilon >0$, we shall explore the methods of Sect. 2 to derive a suitable polynomial p so that

$$\begin{aligned} \Vert \varphi (t,u,0) -f\Vert _\infty = \Vert p(A)b -f\Vert _\infty < \varepsilon . \end{aligned}$$

Starting point: Known sufficient conditions

As a starting point, we take the necessary and sufficient conditions for uniform ensemble reachability developed in [14, Theorem 4, Corollary 3]. Since these are essential, we provide a short recap. That is, if the pair (A, b) is uniformly ensemble reachable, the following necessary conditions are satisfied:

• (N1) The pair $(A(\theta ),b(\theta ))$ is reachable for every $\theta \in {\textbf{P}}$.

• (N2) For any pair of distinct parameters $\theta ,\theta '\in {\textbf{P}}$, the spectra of $A(\theta )$ and $A(\theta ')$ are disjoint:

  $$\begin{aligned} \sigma (A(\theta ))\cap \sigma (A(\theta '))=\emptyset . \end{aligned}$$

Further, (A, b) is uniformly ensemble reachable if it satisfies (N1), (N2) and one of the following sufficiency conditions:

• (S1) For each $\theta \in {\textbf{P}}$ the characteristic polynomial of $A(\theta )$ takes the form

  $$\begin{aligned} z^n - (a_{n-1}z^{n-1} + \cdots +a_1z + a_0(\theta )) \end{aligned}$$

  for some $a_{n-1},...,a_1 \in {\mathbb {C}}$ and $a_0 \in C({\textbf{P}})$.

• (S2) $A(\theta )$ has simple eigenvalues for each $\theta \in {\textbf {P}}$.

A few comments are in order. First, we note that the conditions are valid for continuous-time systems and discrete-time systems. Also, we note that a pointwise characterization for uniform ensemble reachability is currently not available. In [14, Theorem 3 (d)] it is shown that for single-input pairs (A, b) a necessary condition for uniform ensemble reachability is that the eigenvalues are simple on an open dense subset of the parameter space ${\textbf{P}}$. In this sense condition (S2) is not far from being necessary. Furthermore, in certain cases, like the controlled harmonic oscillator, cf. [28], it happens that sufficiency conditions (S1) and (S2) are satisfied at the same time. Besides, to the best of the author’s knowledge, other general sufficient conditions are not known. Moreover, we note that, due to the condition (N2) the function $a_0$ in (S2) is necessarily injective. Thus, $a_0:{\textbf{P}}\rightarrow a_0({\textbf{P}})$ is one-to-one and onto and so $a_0({\textbf{P}})$ defines a Jordan arc. It will turn out that in the special case, where $a_0({\textbf{P}})$ is a subset of the real axis, based on Theorem 1, a constructive method can be obtained by using Bernstein polynomials. The general case is covered by Walsh’s theorem 4. This method turns out to be more involved and will require also Theorems 2, 3, 6 and 7.

Before we begin with the constructions, we note that we start the constructions based on the conditions (N1), (N2) in (S1) or (S2) and discuss additional required conditions along the lines. The corresponding results are formulated subsequently to each construction. This has the advantage that the results and the corresponding polynomials can be spotted simultaneously.

3.1 Methods for condition (S1)

In this subsection, we will explore the circumstances such that the verification of the sufficient condition (S1) (cf. proof of Theorem 4 in [14]) gives rise to methods for the construction of suitable $T>0$ and inputs $u=(u_0,...,u_{T-1})$. Let’s we assume that the pair (A, b) satisfies the conditions (N1), (N2) and (S1). Also, let $f \in {\text {C}}_n({\textbf{P}})$ and $\varepsilon >0$ be given. By condition (N1), the matrix

$$\begin{aligned} T(\theta ) = \begin{pmatrix} b(\theta )&A(\theta )b(\theta )&\cdots&A(\theta )^{n-1}b(\theta ) \end{pmatrix} \end{aligned}$$

is continuously invertible for every $\theta \in {\textbf{P}}$ and one has

$$\begin{aligned} {\tilde{A}}(\theta )= T(\theta )^{-1}A(\theta )T(\theta ) = \left( \begin{array}{ccccc} 0 &{}\quad \hdots &{}\quad \hdots &{}\quad 0 &{}\quad a_0(\theta )\\ 1 &{}\quad \ddots &{}\quad &{}\quad \vdots &{}\quad a_1\\ 0 &{}\quad \ddots &{}\quad \ddots &{}\quad \vdots &{}\quad \vdots \\ \vdots &{}\quad \ddots &{}\quad \ddots &{}\quad 0 &{}\quad \vdots \\ 0 &{}\quad \hdots &{}\quad 0 &{}\quad 1 &{}\quad a_{n-1} \end{array}\right) ,\qquad {\tilde{b}}= T(\theta )^{-1}b(\theta )= \begin{pmatrix} 1\\ 0\\ \vdots \\ 0 \end{pmatrix}, \end{aligned}$$

where $a_0(\theta ),a_1,...,a_{n-1}$ denote the coefficients of the characteristic polynomial of $A(\theta )$. Writing the solution formula in terms of a polynomial p, we get

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty }&= \Vert p(A)b - f\Vert _{\infty } = \Vert T ( p({\tilde{A}}){\tilde{b}} - T^{-1}f)\Vert _{\infty } \\ {}&\le \Vert T\Vert _{M,\infty } \,\Vert p({\tilde{A}}){\tilde{b}} - T^{-1}f\Vert _{\infty }, \end{aligned}$$

where $\Vert \cdot \Vert _M$ is a matrix norm that is submultiplicative to $\Vert \cdot \Vert$ and

$$\begin{aligned} \Vert T\Vert _{M,\infty }:= \sup _{\theta \in {\textbf{P}}} \Vert T(\theta )\Vert _M. \end{aligned}$$

As shown in [14], for

$$\begin{aligned} p(z)= \sum _{k=1}^n p_k\big (z^n-a_{n-1}z^{n-1}-\cdots - a_1z\big )\, z^{k-1}{,} \end{aligned}$$

(18)

we get

$$\begin{aligned} p({\tilde{A}}(\theta )) {\tilde{b}} =\sum _{k=1}^n p_k \big ( A(\theta )^n -a_{n-1} A(\theta )^{n-1} -\cdots {-}a_1 A(\theta ) \big )A(\theta )^{k-1}{\tilde{b}} =\begin{pmatrix} p_1(a_0(\theta ))\\ \vdots \\ p_n(a_0(\theta )) \end{pmatrix} \end{aligned}$$

and, hence,

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty } \le \Vert T\Vert _{M,\infty } \,\sup _{k=1,...,n} \sup _{\theta \in {\textbf{P}}} | p_k(a_0(\theta )) - (T^{-1}f)_k(\theta )| , \end{aligned}$$

where $(T^{-1}f)_k(\theta )$ denotes the kth component of the vector-valued function $\theta \mapsto T^{-1}(\theta )f(\theta )\in {\mathbb {C}}^n$. Moreover, by (N2) and (S1), the function $\theta \mapsto a_0(\theta )$ is one-to-one. Consequently, we get

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty } < \varepsilon {,} \end{aligned}$$

if we can find polynomials $p_1,....p_n$ satisfying

$$\begin{aligned} \sup _{z \in a_0({\textbf{P}})} | p_k(z) - (T^{-1}f)_k \circ a_0^{-1}(z)| < \frac{\varepsilon }{\Vert T\Vert _{M,\infty }} \quad \text { for all } \, k=1,...,n. \end{aligned}$$

(19)

Depending on the structure of $a_0({\textbf{P}})$, we will present two constructive methods.

3.1.1 Case 1: $a_0({\textbf{P}})=[a,b]\subset {\mathbb {R}}$ defines a compact interval

We assume that for each $k=1,...,n$ the function $z \mapsto (T^{-1}f)_k\circ a_0^{-1}(z)$ is Lipschitz continuous³ and let $L_k{>}0$ denote the Lipschitz constant. By Theorem 1, condition (19) is fulfilled if we take $p_k$ as the Bernstein polynomial, cf. (4), corresponding to the function $z \mapsto (T^{-1}f)_k\circ a_0^{-1}(z)$ with degree $n_k \ge 3$ satisfying

$$\begin{aligned} \sqrt{2} \left( 4\,M_{(T^{-1}f)_k\circ a_0^{-1}} + \frac{(b-a)\,L_{k}}{2} \right) {\sqrt{\frac{\log n_k}{n_k}}} < \frac{\varepsilon }{\Vert T\Vert _{M,\infty }}. \end{aligned}$$

Denoting this Bernstein polynomial by $B_{(T^{-1}f)_k\circ a_0^{-1},n_k}$, we conclude that for

$$\begin{aligned} p(z)= \sum _{k=1}^n B_{(T^{-1}f)_k\circ a_0^{-1},n_k}\big (z^n-a_{n-1}z^{n-1}-\cdots - a_1z\big )\, z^{k-1} \end{aligned}$$

it holds

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty } \le \Vert T\Vert _{M,\infty }\, \Vert p({\tilde{A}}) {\tilde{b}}- T^{-1}f\Vert _\infty < \varepsilon . \end{aligned}$$

Considering the monomial representation

$$\begin{aligned} p(z)= \sum _{k=1}^n B_{(T^{-1}f)_k\circ a_0^{-1},n_k}\big (z^n-a_{n-1}z^{n-1}-\cdots - a_1z\big )\, z^{k-1} =\sum _{l=0}^N \alpha _l z^l{,} \end{aligned}$$

(20)

we get the following result.

Theorem 8

Let $(A,b) \in C_{n,n}({\textbf{P}})\times C_n({\textbf{P}})$ satisfy (N1), (N2) and (S1). Let $\varepsilon >0$ and suppose that $f \in C_n({\textbf{P}})$ is such that $(T^{-1}f)_k \circ a_0^{-1} \in {\text {Lip}}(a_0({\textbf{P}}))$ for all $k=1,...,n$ and assume that $a_0({\textbf{P}})$ is a compact real interval. Then, for

$$\begin{aligned} T=N+1 \qquad \text {and} \qquad u=(u_0,...,u_{T-1}) = (\alpha _N,...,\alpha _{0}) \end{aligned}$$

defined in (20), one has

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty } < \varepsilon . \end{aligned}$$

In particular, the result is true if $(A,b) \in C^1_{n,n}({\textbf{P}})\times C^1_n({\textbf{P}})$ and $f\in {\text {Lip}}_n({\textbf{P}})$.

Next we consider the more general problem of constructing polynomials $p_k$ satisfying (19) if $a_0({\textbf{P}})$ defines a Jordan arc. It turns out that this is more involved. We shall see how the proof of Theorem 4 gives rise to a constructive method.

3.1.2 Case 2: $a_0({\textbf{P}})$ defines a Jordan arc

Without loss of generality⁴ we assume that $a_0({\textbf{P}})$ does not contain the origin. In a first step we extent $a_0({\textbf{P}})$ to a Jordan curve $\varGamma _0$ such that the origin lies inside the Jordan curve. Then, take a continuous mapping $\varPhi _0: \overline{ {\text {int}}\varGamma _0} \rightarrow \overline{{\mathbb {D}}}$ which is conformal on ${\text {int}}\varGamma _0$ so that $\varPhi _0(0)=0$. The existence of such a mapping is ensured by Caratheodory’s Theorem, cf. [26, Thm. 13.2.3].

Step 1. Extension of the target function to the Jordan curve

Here, we extend the functions $z \mapsto (T^{-1}f)_k\circ a_0^{-1}(z)$ to $\varGamma _0$, where $T\in C_{n,n}({\textbf{P}})$ denotes the changes of coordinates introduced at the beginning of Sect. 3.1. This can be done as follows. For simplicity, we do not introduce a new notation. For all $z \in \varGamma _0 {\setminus } a_0({\textbf{P}})$ we define

$$\begin{aligned} (T^{-1} f)_k \circ a_0^{-1}(z):= w_{k,1} + (w_{k,2}-w_{k,1})\frac{z-z_{k,1}}{z_{k,2}-z_{k,1}}, \end{aligned}$$

where $w_{k,1}$ and $w_{k,2}$ are the values of $(T^{-1}f)_k\circ a^{-1}_0$ at the end-points $z_{k,1}$ and $z_{k,2}$ of $a_0({\textbf{P}})$, respectively. As a result, we obtain a continuous function

$$\begin{aligned} h_k:= (T^{-1}f)_k\circ a_0^{-1} \circ \varPhi _0^{-1} :\partial {\mathbb {D}} \rightarrow {\mathbb {C}}. \end{aligned}$$

In the following, we assume that $a_0 \in C^2({\textbf{P}})$ and ${\textbf{P}}$ has a parametrization which is twice continuously differentiable. Also we suppose that the extension $\varGamma _0$ has the same properties. Then, by Kellogg’s theorem [25, Theorem 4.3], the functions $h_k$ satisfy a Lipschitz condition provided $(T^{-1}f)_k$ is Lipschitz continuous⁵. We assume that this is the case and let $L_k>0$ denote the corresponding Lipschitz constant.

Step 2. Ansatz for the polynomial p

To derive a suitable polynomial p, we take again (cf. (18))

$$\begin{aligned} p(z)= \sum _{k=1}^n p_k\big (z^n-a_{n-1}z^{n-1}-\cdots - a_1z\big )\, z^{k-1}. \end{aligned}$$

Based on the exploration in Sect. 2.3 (in particular (17)), we will construct the polynomials $p_1,...,p_n$ as follows. For each $k =1,...,n$ we consider the following Fejér-like polynomial

$$\begin{aligned} p_k(z) = \sum _{l=0}^{n_k-1} c_l \, \big (p_{k,1}(z)\big )^l + \sum _{l=1}^{n_k-1} c_{-l} \, \big (p_{k,2}(z) + p_{k,3}(z)\big )^l, \end{aligned}$$

where

$$\begin{aligned} c_l:= \frac{n_k - |l|}{n_k}\, \widehat{h_k} (l),\qquad \, l = -n_k+1,....,n_k-1 \end{aligned}$$

denote the Fourier coefficients of $h_k$, cf. (5). We proceed with a simple but instructive step. By adding and subtracting appropriate terms and using the triangle inequality, one has

$$\begin{aligned} | (T^{-1}f)_k\circ a_0^{-1}(z)&- p_k(z)| \\&\le \left| (T^{-1}f)_k\circ a_0^{-1} (z) - \sum _{l=0}^{n_k-1} c_l \, \big (\varPhi _0(z)\big )^l - \sum _{l=1}^{n_k-1} c_{-l} \, \big (\tfrac{1}{\varPhi _0(z)}\big )^l\right| \\&\quad + \left| \sum _{l=0}^{n_k-1} c_l \, \big (\varPhi _0(z)^l- p_{k,1}(z)^l\big ) \right| \\ {}&\quad + \left| \sum _{l=1}^{n_k-1} c_{-l} \, \Big ( \, \big ( \tfrac{a_{-1}}{z} +p_{k,2}(z)\big )^l - \big (\tfrac{1}{\varPhi _0(z)}\big )^l \Big ) \right| \\&\quad + \left| \sum _{l=1}^{n_k-1} c_{-l} \, \Big ( \, \big (p_{k,2}(z) + p_{k,3}(z) \big )^l - \big ( \tfrac{a_{-1}}{z} +p_{k,2}(z)\big )^l \Big ) \right| , \end{aligned}$$

where

$$\begin{aligned} a_{-1} = \frac{1}{2\pi i} \int _{\partial B_r(0)} \frac{1}{\varPhi _0(\xi )} {\text {d}}\!\xi , \qquad r >0 \, \, \text {such that } \overline{B_r(0)} \subset {\text {int}}\varGamma _0 \end{aligned}$$

denotes the residue of $\tfrac{1}{\varPhi _0}$ at zero.

Step 3. Derivation of $n_k \in {\mathbb {N}}$

Applying Weierstrass second Theorem 2 (b) to the function

$$\begin{aligned} h_k = (T^{-1}f)_k\circ a_0^{-1} \circ \varPhi _0^{-1} :\partial {\mathbb {D}} \rightarrow {\mathbb {C}}{,} \end{aligned}$$

we get

$$\begin{aligned} \sup _{w \in \partial {\mathbb {D}}}\left| h_k (w) - \sum _{l=0}^{n_k-1} c_l \, w^l - \sum _{l=1}^{n_k-1} c_{-l} \, \big (\tfrac{1}{w}\big )^l\right| < 2\sqrt{2}\pi L_{k} \frac{\ln n_k }{n_k} {.} \end{aligned}$$

As $\varPhi _0$ is one-to-one and onto and taking $n_k \in {\mathbb {N}}$ so that

$$\begin{aligned} 2\sqrt{2}\pi L_{k} \frac{\ln n_k }{n_k} \le \frac{\varepsilon }{4\,\Vert T\Vert _{M,\infty }}, \end{aligned}$$

we have

$$\begin{aligned} \sup _{z \in \varGamma _0} \left| (T^{-1}f)_k\circ a_0^{-1} (z) - \sum _{l=0}^{n_k-1} c_l \, \big (\varPhi _0(z)\big )^l - \sum _{l=1}^{n_k-1} c_{-l} \, \big (\tfrac{1}{\varPhi _0(z)}\big )^l\right| < \frac{\varepsilon }{4 \, \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Hence it remains to determine the polynomials $p_{k,1},p_{k,2},p_{k,3}$. This will be carried out in the following three steps. To this end, we pick $\delta _{k}>0$ so that

$$\begin{aligned} \sum _{l=0}^{n_k-1} \left( (1+\delta _k)^l -1 \right) < \frac{\varepsilon }{4\,c_{\max } \Vert T\Vert _{M,\infty }}, \qquad c_{\max } := \max \{{\big |}c_{-n_k+1}{\big |},...,{\big |}c_{n_k-1}{\big |}\}. \end{aligned}$$

Step 4. Derivation of $p_{k,1}$

According to Courant’s Theorem 5, we take a Jordan curve $\varGamma _1$ such that $\varGamma _1 \subset {\mathbb {C}}{\setminus } \overline{ {\text {int}}\varGamma _0}$ and we choose a continuous mapping $\varPhi _1: \overline{ {\text {int}}\varGamma _1} \rightarrow {\mathbb {D}}$ which is conformal on ${\text {int}}\varGamma _1$ satisfying $\varPhi _1 (0)=0$ and

$$\begin{aligned} | \varPhi _0 (z)- \varPhi _1 (z)| < \frac{\delta _k}{2} \qquad \text { for all } z \in \overline{ {\text {int}}\varGamma _0}. \end{aligned}$$

(21)

The construction of $\varGamma _1$ and hence of $\varPhi _1$ depends on the properties of $a_0$. We note that numerical methods for conformal mappings are a research area itself and, thus, it is beyond the scope of this paper. We refer to the comments provided at the end of Sect. 2.

To derive the polynomial $p_{k,1}$ we apply Runge’s little Theorem 3 to the holomorphic function $\varPhi _1: {\text {int}}\varGamma _1 \rightarrow {\mathbb {D}}$ so that

$$\begin{aligned} |p_{k,1}(z) - \varPhi _1(z)|< \frac{\delta _k}{2} \qquad \text { for all } z \in \overline{ {\text {int}}\varGamma _0}. \end{aligned}$$

Then, for all $z \in \overline{ {\text {int}}\varGamma _0}$, one has

$$\begin{aligned} |p_{k,1}(z) - \varPhi _0(z)| \le |p_{k,1}(z) - \varPhi _1(z)| + | \varPhi _1(z) - \varPhi _0(z)|< \delta _k . \end{aligned}$$

Further, as shown in the proof of Theorem 6 (cf.(14)), for all $z\in \varGamma _0$ it follows

$$\begin{aligned} \left| \sum _{l=0}^{n_k-1} c_l \, \big (\varPhi _0(z)^l- p_{k,1}(z)^l\big ) \right| <\frac{\varepsilon }{4 \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Step 5. Derivation of $p_{k,2}$

We consider the mapping $\varPhi _0:\overline{ {\text {int}}\varGamma _0} \rightarrow \overline{{\mathbb {D}}}$ and the residue of $\tfrac{1}{\varPhi _0}$ at zero, i.e.,

$$\begin{aligned} a_{-1} = \frac{1}{2\pi i} \int _{\partial B_r(0)} \frac{1}{\varPhi _0(\xi )} {\text {d}}\!\xi , \qquad r >0 \, \, \text {such that } \overline{B_r(0)} \subset {\text {int}}\varGamma _0. \end{aligned}$$

Moreover, let $\Psi _0: {\text {int}}\varGamma _0\rightarrow {\mathbb {C}}$,

$$\begin{aligned} \Psi _0(z) = {\left\{ \begin{array}{ll} \frac{1}{\varPhi _0(z)} - \frac{a_{-1}}{z} &{} z \in {\text {int}}\varGamma _0\setminus \{0\}\\ \frac{1}{2\pi i} \int _{\partial B_r(0)} \frac{1}{\xi \varPhi _0(\xi )} {\text {d}}\!\xi &{} z=0, \end{array}\right. } \end{aligned}$$

(22)

denote the holomorphic extension, which is continuous on $\overline{ {\text {int}}\varGamma _0}$. To determine $p_{k,2}$, we take another Jordan curve $\varGamma _2$ such that $\varGamma _2 \subset {\mathbb {C}}{\setminus } \overline{ {\text {int}}\varGamma _0}$ and we choose a continuous mapping $\varPhi _2: \overline{ {\text {int}}\varGamma _2} \rightarrow {\mathbb {D}}$ which is conformal on ${\text {int}}\varGamma _2$ satisfying $\varPhi _2 (0)=0$ and

$$\begin{aligned} | \varPhi ^{-1}_0 (w) - \varPhi ^{-1}_2 (w)| < \frac{\delta _k}{2L_{\Psi _0}} \qquad \text { for all } w \in \overline{{\mathbb {D}}}, \end{aligned}$$

(23)

where $L_{\Psi _0}>0$ denotes the Lipschitz constant of $\Psi _0$.⁶ Then, for all $z =\varPhi _2^{-1}(w) \in \overline{ {\text {int}}\varGamma _0}$, we have

$$\begin{aligned} | \Psi _0 (z)- \Psi _0\circ \varPhi ^{-1}_0\circ \varPhi _2 (z)| \le L_{\Psi _0} \, | \varPhi ^{-1}_2 (w) - \varPhi ^{-1}_0 (w)| < \frac{\delta _k}{2}. \end{aligned}$$

Next, we apply Runge’s little theorem to the holomorphic function

$$\begin{aligned} \Psi _0\circ \varPhi ^{-1}_0\circ \varPhi _2: {\text {int}}\varGamma _2 \rightarrow {\mathbb {C}}\end{aligned}$$

and determine the polynomial $p_{k,2}$ such that

$$\begin{aligned} | p_{k,2} (z) - \Psi _0\circ \varPhi ^{-1}_0\circ \varPhi _2 (z)| < \frac{\delta _k}{2} \qquad \text { for all } z \in \overline{ {\text {int}}\varGamma _0}, \end{aligned}$$

Thus, as shown in the proof of Theorem 6 (similar to (14)), for all $z\in \varGamma _0$ it follows

$$\begin{aligned} \left| \sum _{l=1}^{n_k-1} c_{-l} \, \Big ( \, \big ( \tfrac{a_{-1}}{z} +p_{k,2}(z)\big )^l - \big (\tfrac{1}{\varPhi _0(z)}\big )^l \Big ) \right| <\frac{\varepsilon }{4 \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Step 6. Derivation of $p_{k,3}$

In this step we will apply Runge’s little theorem to the holomorphic function $z \mapsto \tfrac{a_{-1}}{z}$ on $a_0({\textbf{P}})$. To this end, let

$$\begin{aligned} \nu := \max _{z \in a_0({\textbf{P}})} \left| \tfrac{1}{z} \right| , \qquad \mu _k:= \max _{z \in a_0({\textbf{P}})} \left| p_{k,2}(z) \right| . \end{aligned}$$

Then, we take $\eta _k>0$ such that

$$\begin{aligned} \eta _k \cdot \frac{1-(2\eta _k + |a_{-1}|\nu +\mu _k)^{n_k-1}}{1-(2\eta _k + |a_{-1}|\nu +\mu _k)} < \frac{\varepsilon }{4 (n_k-1)c_{\max } \Vert T\Vert _{M,\infty }}. \end{aligned}$$

The polynomial $p_{k,3}$ is then obtained by applying Runge’s little theorem such that for all $z \in a_0({\textbf{P}})$ one has

$$\begin{aligned} |\tfrac{a_{-1}}{z} -p_{k,3}(z)|< \eta _k . \end{aligned}$$

Then, as shown in the proof of Theorem 4, we get

$$\begin{aligned} \left| \sum _{l=1}^{n_k-1} c_{-l} \, \Big ( \, \big (p_{k,2}(z) + p_{k,3}(z) \big )^l - \big ( \tfrac{a_{-1}}{z} +p_{k,2}(z)\big )^l \Big ) \right| < \frac{\varepsilon }{4 \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Thus, we have

$$\begin{aligned} | (T^{-1}f)_k\circ a_0^{-1}(z) - p_k(z)| < \frac{\varepsilon }{ \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Step 7. Monomial representation of p

Finally, we determine the monomial representation of the polynomial

$$\begin{aligned} p(z) = \sum _{k=1}^n p_k(z^n+a_{n-1}z^{n-1}+ \cdots + a_1\,z) z^{k-1} =:\sum _{l=0}^N \beta _l z^l. \end{aligned}$$

(24)

This construction yields the following

Theorem 9

Let ${\textbf{P}}$ be a Jordan arc having a twice continuously differentiable parametrization. Assume that $(A,b)\in C^1_{n,n}({\textbf{P}})\times C_n^1({\textbf{P}})$ satisfies (N1), (N2) and (S1). Let $\varepsilon >0$, $f \in {\text {C}}_n({\textbf{P}})$ and $a_0 \in C^2({\textbf{P}})$ be a Jordan arc in the complex plane such that $(T^{-1}f)_k \circ a_0^{-1} \in {\text {Lip}}(a_0({\textbf{P}}))$ for all $k=1,...,n$. Then, for

$$\begin{aligned} T=N+1 \qquad \text {and} \qquad u=(u_0,...,u_{T-1}) = (\beta _N,...,\beta _{0}) \end{aligned}$$

defined in (24), one has

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty } < \varepsilon . \end{aligned}$$

In particular, the result is true if $(A,b) \in C^1_{n,n}({\textbf{P}})\times C^1_n({\textbf{P}})$ and $f\in {\text {Lip}}_n({\textbf{P}})$.

3.2 Methods for condition (S2)

In this section we consider methods to construct suitable inputs based on the necessary conditions (N1), (N2) and the sufficiency condition (S2). Let $f \in {\text {C}}_n({\textbf{P}})$ and $\varepsilon >0$ be given. After applying a change of coordinates $T(\theta )$, we consider the pair

$$\begin{aligned} {\tilde{A}}(\theta ):= T(\theta )^{-1} A(\theta )T(\theta )=\begin{pmatrix} \lambda _1(\theta ) &{} &{} \\ &{}\quad \ddots &{}\quad \\ &{}\quad &{}\quad \lambda _n(\theta ) \end{pmatrix},\quad {\tilde{b}}(\theta ):=T(\theta )^{-1}b(\theta ) = \begin{pmatrix} 1\\ \vdots \\ 1 \end{pmatrix}, \end{aligned}$$

where $\lambda _1,...,\lambda _n$ denote the distinct eigenvalue Jordan arcs. Writing the solution formula in terms of a polynomial p, we get

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty }&= \Vert p(A)b - f\Vert _{\infty } = \Vert T ( p({\tilde{A}}){\tilde{b}} - T^{-1}f)\Vert _{\infty } \\ {}&\le \Vert T\Vert _{M,\infty } \,\Vert p({\tilde{A}}){\tilde{b}} - T^{-1}f\Vert _{\infty }\\&= \Vert T\Vert _{M,\infty } \,\left\| \begin{pmatrix} p(\lambda _1) - (T^{-1}f)_1\\ \vdots \\ p(\lambda _n) - (T^{-1}f)_n \end{pmatrix} \right\| _{\infty }, \end{aligned}$$

where $(T^{-1}f)_k$ denotes the kth component of the vector-valued function $\theta \mapsto T^{-1}(\theta )f(\theta )\in {\mathbb {C}}^n$. Thus, we conclude that

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty } < \varepsilon {,} \end{aligned}$$

if we can find a polynomial p so that for all $k=1,...,n$ it holds

$$\begin{aligned} \sup _{z \in \lambda _k({\textbf{P}})} | p(z) - (T^{-1}f)_k \circ \lambda _k^{-1}(z)| < \frac{\varepsilon }{\Vert T\Vert _{M,\infty }}. \end{aligned}$$

Compared to the methods for (S1), we are a looking for a single polynomial working for every component of the target function. With other words, the polynomial p solves n approximation problems simultaneously, where each problem is defined on an individual domain. Inspired by [3], using the ansatz

$$\begin{aligned} p(z) = \sum _{k=1}^n p_k(z)q_k(z), \end{aligned}$$

we will split the problem into n subproblems. Postponing the technicalities for a moment, on the one hand we will determine $p_1,...,p_n$ so that for each $k\in \{1,...,n\}$ the polynomial $p_k$ approximates the function $(T^{-1}f)_k \circ \lambda _k^{-1}$ over the image of $\lambda _k:{\textbf{P}}\rightarrow {\mathbb {C}}$. Let $\lambda _k({\textbf{P}})$ denote the image of $\lambda _k$. On the other hand, the polynomials $q_1,...,q_n$ are derived by approximating the indicator functions

$$\begin{aligned} h_k:\lambda ({\textbf{P}}) \rightarrow {\mathbb {C}}, \quad h_k(z)= {\left\{ \begin{array}{ll} 1 &{} \text { if } z \in \lambda _k ({\textbf{P}}) \\ 0 &{} \text { if } z \in \lambda ({\textbf{P}}) \setminus \lambda _k({\textbf{P}}), \end{array}\right. } \end{aligned}$$

where

$$\begin{aligned} \lambda ({\textbf{P}}) = \bigcup _{k=1}^n \lambda _k({\textbf{P}}). \end{aligned}$$

Note that the assumptions imply that the sets $\lambda _1({\textbf{P}}),..., \lambda _n({\textbf{P}})$ are pairwise disjoint and so $h_1,...,h_n$ are holomorphic functions.

A precise outline for the construction of the polynomials $p_1,...,p_n$ and $q_1,...,q_n$ is presented in the following. In doing so, we will distinguish two cases. Exemplary we will discuss the construction process for polynomials $p_k$ and $q_k$, for some $k \in \{1,...,n\}$.

3.2.1 Case 1: $\lambda _k({\textbf{P}})=[a,b]\subset {\mathbb {R}}$ defines a compact interval

We suppose the function $z \mapsto (T^{-1}f)_k\circ \lambda _k^{-1}(z)$ is Lipschitz continuous.⁷. Let $L_k{>}0$ denote the Lipschitz constant. Then, by Theorem 1 the condition

$$\begin{aligned} \sup _{z \in \lambda _k({\textbf{P}})} | p_k(z) - (T^{-1}f)_k \circ \lambda _k^{-1}(z)| < \frac{\varepsilon }{\Vert T\Vert _{M,\infty }} \end{aligned}$$

(25)

is achieved by taking $p_k= B_{(T^{-1}f)_k\circ \lambda _k^{-1},n_k}$ to be the Bernstein polynomial of the function $z \mapsto (T^{-1}f)_k\circ \lambda _k^{-1}(z)$ of degree $n_k \ge 3$ satisfying

$$\begin{aligned} \sqrt{2} \left( 4\,M_{(T^{-1}f)_k\circ \lambda _k^{-1}} + \frac{(b-a)\,L_{k}}{2} \right) {\sqrt{\frac{\log n_k}{n_k}}} < \frac{\varepsilon }{\Vert T\Vert _{M,\infty }}. \end{aligned}$$

After we determined $p_1,...,p_n$, we define

$$\begin{aligned} \alpha _{k,l}:= \sup _{\theta \in {\textbf{P}}}|B_{(T^{-1}f)_l\circ \lambda _l^{-1},n_l}(\lambda _k(\theta ))|. \end{aligned}$$

The polynomial $q_k$ is derived by applying Runge’s little theorem to the holomorphic function $h_k$ such that for each $k=1,...,n$ it holds

$$\begin{aligned} \sup _{z \in \lambda _k({\textbf{P}})} | q_k(z) - h_k(z))| < \frac{\varepsilon }{3\Vert T\Vert _{M,\infty } \, \sum _{l=1}^n \alpha _{k,l}}. \end{aligned}$$

(26)

For every $k=1,...,n$, based on conditions (25) and (26), we get

$$\begin{aligned} | (T^{-1}f)_k \circ \lambda _k^{-1}(z) -p(z)|&\le | (T^{-1}f)_k \circ \lambda _k^{-1}(z)- p_k(z)q_k(z)| + \sum _{l\ne k} |q_l(z)|\, |p_l(z)|\\&\le | (T^{-1}f)_k \circ \lambda _k^{-1}(z)- p_k(z)| + |p_k(z)||1-q_k(z)| \\&\quad + \sum _{l\ne k} |q_l(z)|\, |p_l(z)|\\&\le \frac{\varepsilon }{3 \Vert T\Vert _{M,\infty }} + \frac{\varepsilon \, \alpha _{k,k}}{3 \Vert T\Vert _{M,\infty }\sum _{l=1}^n \alpha _{k,l}}\\&\quad + \sum _{l \ne k} \frac{\varepsilon \, \alpha _{k,l}}{3 \Vert T\Vert _{M,\infty } \sum _{l=1}^n \alpha _{k,l}} \le \frac{\varepsilon }{ \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Consequently, by considering the monomial representation

$$\begin{aligned} p(z)= \sum _{k=1}^n q_k(z)\, B_{(T^{-1}f)_k\circ \lambda _k^{-1},n_k} (z) =\sum _{k=0}^N \beta _k z^k \end{aligned}$$

(27)

we get the following result.

Theorem 10

Let $(A,b) \in C_{n,n}({\textbf{P}})\times C_n({\textbf{P}})$ satisfy (N1), (N2) and (S2). Let $\varepsilon >0$ and suppose that $f \in C_n({\textbf{P}})$ is such that $(T^{-1}f)_k \circ \lambda _k^{-1} \in {\text {Lip}}(\lambda _k({\textbf{P}}))$ for all $k=1,...,n$ and assume that $\lambda _1({\textbf{P}}),...,\lambda _n({\textbf{P}})$ are compact real intervals. Then, for

$$\begin{aligned} T=N+1 \qquad \text {and} \qquad u=(u_0,...,u_{T-1}) = (\beta _N,...,\beta _{0}) \end{aligned}$$

defined in (27), one has

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty } < \varepsilon . \end{aligned}$$

In particular, the result is true if $(A,b) \in C^1_{n,n}({\textbf{P}})\times C^1_n({\textbf{P}})$ and $f\in {\text {Lip}}_n({\textbf{P}})$.

Next we consider the more general case where the eigenvalue functions define Jordan arcs.

3.2.2 Case 2: $\lambda _1({\textbf{P}}),...,\lambda _n({\textbf{P}})$ define Jordan arcs

The construction process is similar to the methods for (S1). Hence, in order to avoid a lengthy presentation we will go through the analogous construction steps more quickly.

Step 1. Extension of the Jordan arcs and the target function

Without loss of generality we will assume that none of the Jordan arcs contains the origin. In a first step we extent the Jordan arcs $\lambda _1({\textbf{P}}),...,\lambda _n({\textbf{P}})$ to Jordan curves $\varGamma _1,..., \varGamma _n$ such that the origin lies inside of each Jordan curve. Then we take a continuous bijective mappings $\varPhi _k: \overline{ {\text {int}}\varGamma _k} \rightarrow \overline{{\mathbb {D}}}$ which are conformal on ${\text {int}}\varGamma _k$ so that $\varPhi _k(0)=0$. The existence of such mappings is ensured by Caratheodory’s Theorem, cf. [26, Thm. 13.2.3].

Analogous to Step 1 in the methods for (S1), we extend the functions $z \mapsto (T^{-1}f)_k\circ \lambda _k^{-1}(z)$ to $\varGamma _k$, where $T\in C_{n,n}({\textbf{P}})$ denotes the changes of coordinates introduced at the beginning of Sect. 3.2. As a result, we obtain continuous functions

$$\begin{aligned} h_k:= (T^{-1}f)_k\circ \lambda _k^{-1} \circ \varPhi _k^{-1} :\partial {\mathbb {D}} \rightarrow {\mathbb {C}}{, \quad k=1,...,n}. \end{aligned}$$

If $\lambda _1,...,\lambda _n \in C^2({\textbf{P}})$ and ${\textbf{P}}$ has a parametrization which is twice continuously differentiable and the extensions $\varGamma _1,...,\varGamma _n$ have the same properties, then by Kellogg’s theorem [25, Theorem 4.3], the functions $h_k$ satisfy a Lipschitz condition provided the $(T^{-1}f)_k$ are Lipschitz continuous⁸. Let $L_k>0$ denote the corresponding Lipschitz constant.

Step 2. Ansatz for the polynomial p

As in the previous case, to derive a suitable polynomial, we set

$$\begin{aligned} p(z)= \sum _{k=1}^n p_k(z)\, q_{k}(z). \end{aligned}$$

Based on the exploration in Sect. 2.3 (in particular (17)), for each $k =1,...,n$ we make the ansatz

$$\begin{aligned} p_k(z) = \sum _{k=0}^{n_k-1} c_l \, \big (p_{k,1}(z)\big )^l + \sum _{l=1}^{n_k-1} c_{-l} \, \big (p_{k,2}(z) + p_{k,3}(z)\big )^l, \end{aligned}$$

where

$$\begin{aligned} c_l:= \frac{n_k - |l|}{n_k}\, \widehat{h_k} (l),\qquad \, l = -n_k+1,....,n_k-1 \end{aligned}$$

denote the Fourier coefficients corresponding to the function $h_k$. Again we take a simple but instructive step

$$\begin{aligned} | (T^{-1}f)_k\circ \lambda _k^{-1}(z)&- p_k(z)| \\&\le \left| (T^{-1}f)_k\circ a_0^{-1} (z) - \sum _{l=0}^{n_k-1} c_l \, \big (\varPhi _k(z)\big )^l - \sum _{l=1}^{n_k-1} c_{-l} \, \big (\tfrac{1}{\varPhi _k(z)}\big )^l\right| \\&\quad + \left| \sum _{l=0}^{n_k-1} c_l \, \big (\varPhi _k(z)^l- p_{k,1}(z)^l\big ) \right| \\&\quad + \left| \sum _{l=1}^{n_k-1} c_{-l} \, \Big ( \, \big ( \tfrac{a_{-1,k}}{z} +p_{k,2}(z)\big )^l - \big (\tfrac{1}{\varPhi _k(z)}\big )^l \Big ) \right| \\ {}&\quad + \left| \sum _{l=1}^{n_k-1} c_{-l} \, \Big ( \, \big (p_{k,2}(z) + p_{k,3}(z) \big )^l - \big ( \tfrac{a_{-1,k}}{z} +p_{k,2}(z)\big )^l \Big ) \right| , \end{aligned}$$

where

$$\begin{aligned} a_{-1,k} = \frac{1}{2\pi i} \int _{\partial B_r(0)} \frac{1}{\varPhi _k(\xi )} {\text {d}}\!\xi , \qquad r >0 \, \, \text {such that } \overline{B_r(0)} \subset {\text {int}}\varGamma _0 \end{aligned}$$

denotes the residue of $\tfrac{1}{\varPhi _k}$ at zero.

Step 3. Derivation of $n_k \in {\mathbb {N}}$

Applying Weierstrass second Theorem 2 (b) to the function

$$\begin{aligned} h_k = (T^{-1}f)_k\circ \lambda _k^{-1} \circ \varPhi _k^{-1} :\partial {\mathbb {D}} \rightarrow {\mathbb {C}} \end{aligned}$$

we get

$$\begin{aligned} \sup _{w \in \partial {\mathbb {D}}}\left| h_k (w) - \sum _{l=0}^{n_k-1} c_l \, w^l - \sum _{l=1}^{n_k-1} c_{-l} \, \big (\tfrac{1}{w}\big )^l\right| < 2\sqrt{2}\pi L_{k} \frac{\ln n_k }{n_k} . \end{aligned}$$

As $\varPhi _k$ is one-to-one and onto and taking $n_k \in {\mathbb {N}}$ so that

$$\begin{aligned} 2\sqrt{2}\pi L_{k} \frac{\ln n_k }{n_k} \le \frac{\varepsilon }{12\,\Vert T\Vert _{M,\infty }}, \end{aligned}$$

we have

$$\begin{aligned} \sup _{z \in \varGamma _k} \left| (T^{-1}f)_k\circ \lambda _k^{-1} (z) - \sum _{l=0}^{n_k-1} c_l \, \big (\varPhi _k(z)\big )^l - \sum _{l=1}^{n_k-1} c_{-l} \, \big (\tfrac{1}{\varPhi _k(z)}\big )^l\right| < \frac{\varepsilon }{12 \, \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Hence, it remains to determine the polynomials $p_{k,1},p_{k,2},p_{k,3}$. This will be carried out in the following three steps. To this end, we pick $\delta _{k}>0$ so that

$$\begin{aligned} \sum _{l=0}^{n_k-1} \left( (1+\delta _k)^l -1 \right) < \frac{\varepsilon }{12\,c_{\max } \Vert T\Vert _{M,\infty }}, \qquad c_{\max } := \max \{{|}c_{-n_k+1}{|},...,{|}c_{n_k-1}{|}\}. \end{aligned}$$

Step 4. Derivation of $p_{k,1}$

According to Courant’s Theorem 5, we take a Jordan curve $\varGamma _{k,1}$ such that $\varGamma _{k,1} \subset {\mathbb {C}}{\setminus } \overline{ {\text {int}}\varGamma _k}$ and we choose a continuous mapping $\varPhi _{k,1}: \overline{ {\text {int}}\varGamma _{k,1}} \rightarrow \overline{{\mathbb {D}}}$ which is conformal on ${\text {int}}\varGamma _{k,1}$ satisfying $\varPhi _{k,1} (0)=0$ and

$$\begin{aligned} | \varPhi _k (z)- \varPhi _{k,1} (z)| < \frac{\delta _k}{2} \qquad \text { for all } z \in \overline{ {\text {int}}\varGamma _k}. \end{aligned}$$

(28)

The construction of $\varGamma _{k,1}$ and hence of $\varPhi _{k,1}$ depends on the properties of $\lambda _k$. Recall that constructive methods for conformal mappings are beyond the scope of this paper. We refer to the comments provided at the end of Sect. 2.

To derive the polynomial $p_{k,1}$, we apply Runge’s little Theorem 3 to the conformal mapping $\varPhi _{k,1}$ so that

$$\begin{aligned} |p_{k,1}(z) - \varPhi _{k,1}(z)|< \frac{\delta _k}{2} \qquad \text { for all } z \in \overline{ {\text {int}}\varGamma _k}. \end{aligned}$$

Then, for all $z \in \overline{ {\text {int}}\varGamma _k}$ one has

$$\begin{aligned} |p_{k,1}(z) - \varPhi _k(z)| \le |p_{k,1}(z) - \varPhi _{k,1}(z)| + | \varPhi _{k,1}(z) - \varPhi _k(z)|< \delta _k . \end{aligned}$$

Further, as shown in the proof of Theorem 6 (cf.(14)), for all $z\in \varGamma _k$ it follows

$$\begin{aligned} \left| \sum _{l=0}^{n_k-1} c_l \, \big (\varPhi _k(z)^l- p_{k,1}(z)^l\big ) \right| <\frac{\varepsilon }{12 \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Step 5. Derivation of $p_{k,2}$

We consider the holomorphic function $\Psi _k: {\text {int}}\varGamma _k\rightarrow {\mathbb {C}}$,

$$\begin{aligned} \Psi _k(z) = {\left\{ \begin{array}{ll} \frac{1}{\varPhi _k(z)} - \frac{a_{-1,k}}{z} &{} z \in {\text {int}}\varGamma _k\setminus \{0\}\\ \frac{1}{2\pi i} \int _{\partial B_r(0)} \frac{1}{\xi \varPhi _k(\xi )} {\text {d}}\!\xi &{} z=0, \end{array}\right. } \end{aligned}$$

(29)

which is continuous on $\overline{ {\text {int}}\varGamma _k}$. To determine $p_{k,2}$, we take another Jordan curve $\varGamma _{k,2}$ such that $\varGamma _{k,2} \subset {\mathbb {C}}\setminus \overline{ {\text {int}}\varGamma _k}$ and we choose a continuous mapping $\varPhi _{k,2}: \overline{ {\text {int}}\varGamma _{k,2}} \rightarrow \overline{{\mathbb {D}}}$ which is conformal on ${\text {int}}\varGamma _{k,2}$ satisfying $\varPhi _{k,2} (0)=0$ and

$$\begin{aligned} | \varPhi ^{-1}_k (w) - \varPhi ^{-1}_{k,2} (w)| < \frac{\delta _k}{2L_{\Psi _k}} \qquad \text { for all } w \in \overline{{\mathbb {D}}}, \end{aligned}$$

(30)

where $L_{\Psi _k}>0$ denotes the Lipschitz constant of $\Psi _k$. Then, for all $z =\varPhi _{k,2}^{-1}(w) \in \overline{ {\text {int}}\varGamma _k}$ we have

$$\begin{aligned} | \Psi _k (z)- \Psi _k\circ \varPhi ^{-1}_0\circ \varPhi _{k,2} (z)| \le L_{\Psi } \, | \varPhi ^{-1}_{k,2} (w) - \varPhi ^{-1}_k (w)| < \frac{\delta _k}{2}. \end{aligned}$$

Next, we apply Runge’s little theorem to the holomorphic function

$$\begin{aligned} \Psi _k\circ \varPhi ^{-1}_k\circ \varPhi _{k,2}: {\text {int}}\varGamma _{k,2} \rightarrow {\mathbb {C}}\end{aligned}$$

and determine the polynomial $p_{k,2}$ such that

$$\begin{aligned} | p_{k,2} (z) - \Psi _k\circ \varPhi ^{-1}_k\circ \varPhi _{k,2} (z)| < \frac{\delta _k}{2} \qquad \text { for all } z \in \overline{ {\text {int}}\varGamma _k}, \end{aligned}$$

Thus, as shown in the proof of Theorem 6 (similar to (14)), for all $z\in \varGamma _k$ it follows

$$\begin{aligned} \left| \sum _{l=1}^{n_k-1} c_{-l} \, \Big ( \, \big ( \tfrac{a_{-1,k}}{z} +p_{k,2}(z)\big )^l - \big (\tfrac{1}{\varPhi _k(z)}\big )^l \Big ) \right| <\frac{\varepsilon }{12 \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Step 6. Derivation of $p_{k,3}$

In this step we will apply Runge’s little theorem to the holomorphic function $z \mapsto \tfrac{a_{-1,k}}{z}$ on $\lambda _k({\textbf{P}})$. More precisely, let

$$\begin{aligned} \nu _k:= \max _{z \in \lambda _k({\textbf{P}})} \left| \tfrac{1}{z} \right| , \qquad \mu _k:= \max _{z \in \lambda _k({\textbf{P}})} \left| p_{k,2}(z) \right| . \end{aligned}$$

Then, we take $\eta _k>0$ such that

$$\begin{aligned} \eta _k \cdot \frac{1-(2\eta _k + |a_{-1,k}|\nu _k +\mu _k)^{n_k-1}}{1-(2\eta _k + |a_{-1,k}|\nu _k +\mu _k)} < \frac{\varepsilon }{12 (n_k-1)c_{\max } \Vert T\Vert _{M,\infty }}. \end{aligned}$$

The polynomial $p_{k,3}$ is then obtained by applying Runge’s little theorem such that for all $z \in \lambda _k({\textbf{P}})$ one has

$$\begin{aligned} |\tfrac{a_{-1,k}}{z} -p_{k,3}(z)|< \eta _k . \end{aligned}$$

Then, as shown in the proof of Theorem 4, we get

$$\begin{aligned} \left| \sum _{l=1}^{n_k-1} c_{-l} \, \Big ( \, \big (p_{k,2}(z) + p_{k,3}(z) \big )^l - \big ( \tfrac{a_{-1}}{z} +p_{k,2}(z)\big )^l \Big ) \right| <\frac{\varepsilon }{12 \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Thus, we have

$$\begin{aligned} | (T^{-1}f)_k\circ \lambda _k^{-1}(z)&- p_k(z)| < \frac{\varepsilon }{3 \Vert T\Vert _{M,\infty }}. \end{aligned}$$

(31)

Step 7. Derivation of $q_{k}$

Let $\lambda _k({\textbf{P}})$ denote the image of $\lambda _k:{\textbf{P}}\rightarrow {\mathbb {C}}$ and

$$\begin{aligned} \lambda ({\textbf{P}}) = \bigcup _{k=1}^n \lambda _k({\textbf{P}}). \end{aligned}$$

The indicator function on $\lambda ({\textbf{P}})$ corresponding to the Jordan arc $\lambda _k$ is denoted by

$$\begin{aligned} h_k:\lambda ({\textbf{P}}) \rightarrow {\mathbb {C}}, \quad h_k(z)= {\left\{ \begin{array}{ll} 1 &{} \text { if } z \in \lambda _k ({\textbf{P}}) \\ 0 &{} \text { if } z \in \lambda ({\textbf{P}}) \setminus \lambda _k({\textbf{P}}). \end{array}\right. } \end{aligned}$$

Recall that the assumptions imply that $\lambda _1({\textbf{P}}),..., \lambda _n({\textbf{P}})$ are pairwise disjoint and so $h_k$ is a holomorphic function. Next, we define

$$\begin{aligned} \alpha _{k,l}:= \sup _{\theta \in {\textbf{P}}}|p_l(\lambda _k(\theta ))|. \end{aligned}$$

The polynomial $q_k$ is then derived by applying Runge’s little theorem to $h_k$ such that

$$\begin{aligned} \sup _{z \in \lambda _k({\textbf{P}})} | q_k(z) - h_k(z))| < \frac{\varepsilon }{3\Vert T\Vert _{M,\infty } \, \sum _{l=1}^n \alpha _{k,l}}. \end{aligned}$$

(32)

Based on conditions (31) and (32), for every $k=1,...,n$ we get

$$\begin{aligned}&| (T^{-1}f)_k \circ \lambda _k^{-1}(z) -p(z)| \\ {}&\quad \le | (T^{-1}f)_k \circ \lambda _k^{-1}(z)- p_k(z)q_k(z)| + \sum _{l\ne k} |q_l(z)|\, |p_l(z)|\\&\quad \le | (T^{-1}f)_k \circ \lambda _k^{-1}(z)- p_k(z)| + |p_k(z)||1-q_k(z)| + \sum _{l\ne k} |q_l(z)|\, |p_l(z)|\\&\quad \le \frac{\varepsilon }{3 \Vert T\Vert _{M,\infty }} + \frac{\varepsilon \, \alpha _{k,k}}{3 \Vert T\Vert _{M,\infty }\sum _{l=1}^n \alpha _{k,l}} + \sum _{l \ne k} \frac{\varepsilon \, \alpha _{k,l}}{3 \Vert T\Vert _{M,\infty } \sum _{l=1}^n \alpha _{k,l}} \le \frac{\varepsilon }{ \Vert T\Vert _{M,\infty }}. \end{aligned}$$

Step 8. Monomial representation of p

Finally, using the polynomials $p_1,...,p_n$ and $q_1,...,q_n$, we determine the monomial representation of the polynomial

$$\begin{aligned} p(z) = \sum _{k=1}^n p_k(z) q_k(z) =:\sum _{l=0}^N \beta _l z^l. \end{aligned}$$

(33)

This construction yields the following result.

Theorem 11

Let ${\textbf{P}}$ be a Jordan arc having a twice continuously differentiable parametrization. Assume that $(A,b)\in C^1_{n,n}({\textbf{P}})\times C_n^1({\textbf{P}})$ satisfies (N1), (N2) and (S1). Let $\varepsilon >0$, $f \in {\text {C}}_n({\textbf{P}})$ and $\lambda _1,...,\lambda _n \in C^2({\textbf{P}})$ be Jordan arcs in the complex plane such that $(T^{-1}f)_k \circ \lambda _k^{-1} \in {\text {Lip}}(\lambda _k({\textbf{P}}))$ for all $k=1,...,n$. Then, for

$$\begin{aligned} T=N+1 \qquad \text {and} \qquad u=(u_0,...,u_{T-1}) = (\beta _N,...,\beta _{0}) \end{aligned}$$

defined in (33), one has

$$\begin{aligned} \Vert \varphi (T,u,0) - f\Vert _{\infty } < \varepsilon . \end{aligned}$$

In particular, the result is true if $(A,b) \in C^1_{n,n}({\textbf{P}})\times C^1_n({\textbf{P}})$ and $f\in {\text {Lip}}_n({\textbf{P}})$.

4 Methods for (S2) for continuous-time single-input systems

In this section, we investigate how to get constructive methods for continuous-time single-input systems. First, we note that a direct application of the discrete-time methods above is not immediately possible. The reason is that, for $u\in L^1([0,T])$ the solution to

$$\begin{aligned} \frac{\partial x}{\partial t} (t,\theta )=A(\theta )x(t,\theta )+b(\theta )u(t), \quad x(0,\theta )=x_0(\theta ) \end{aligned}$$

is given by

$$\begin{aligned} \varphi (T,u,x_0) (\theta ) = \textrm{e}^{TA(\theta )} x_0(\theta ) + \int _0^T \textrm{e}^{(T-s)A(\theta )}b(\theta )u(s){\text {d}}\!s, \end{aligned}$$

which is not of the form $p(A(\theta ))b(\theta )$ for some polynomial p. In the following, we present an approach that enables us to trace the input construction problem for continuous-time systems back to the methods just presented in the previous section. The basic idea simply is to approximate the integral in the solution formula by a polynomial using piecewise constant inputs. The analysis will be carried out in six steps. Before we start, there are two comments in order. For continuous-time systems, the final time can be chosen arbitrarily, i.e., we can fix some $T>0$ in advance. As a consequence, for continuous-time systems we are not restricted to start at the origin. Indeed, let $x_0 \in C_n({\textbf{P}})$ denote the initial states and let $f\in C_n({\textbf{P}})$ denote the final states, then using ${\tilde{f}} \in C_n({\textbf{P}})$, defined by

$$\begin{aligned} {\tilde{f}}(\theta ):= f(\theta ) - \textrm{e}^{TA(\theta )} x_0(\theta ) \end{aligned}$$

an input $u \in L^1([0,T])$ satisfies $\Vert \varphi (T,u,x_0)-f\Vert _{\infty }<\varepsilon$ if and only if $\Vert \varphi (T,u,0)-{\tilde{f}}\Vert _{\infty }<\varepsilon$.

Suppose the pair (A, b) satisfies the conditions (N1), (N2) and (S2) and let $f \in C_n({\textbf{P}})$ and $\varepsilon >0$ be given. As $T>0$ can chosen arbitrarily, we set $T = N\tau$ for some $\tau >0$ and some $N \in {\mathbb {N}}$ that will be specified in the construction process.

Step 1: Diagonalize $A(\theta )$

By assumptions (N2) and (S2) there is a continuous change of coordinates $T(\theta )$ such that

$$\begin{aligned} T(\theta )^{-1} A(\theta ) T(\theta ) = \begin{pmatrix} \lambda _1(\theta ) &{}\quad &{}\quad 0\\ &{}\quad \ddots &{}\quad \\ 0 &{}\quad &{}\quad \lambda _n(\theta ) \end{pmatrix} \quad \text { and } \quad T(\theta )^{-1} b(\theta )= \begin{pmatrix} 1\\ \vdots \\ 1 \end{pmatrix}. \end{aligned}$$

The solution can be written as

$$\begin{aligned} \varphi (N\tau ,u,0) (\theta ) = T(\theta ) \begin{pmatrix} {\tilde{\varphi }}_1(N\tau ,u,0)(\theta ) \\ \vdots \\ {\tilde{\varphi }}_n(N\tau ,u,0)(\theta ) \end{pmatrix}, \end{aligned}$$

where

$$\begin{aligned} {\tilde{\varphi }}_k(N\tau ,u,0)(\theta )= \int _0^{N\tau } \textrm{e}^{(N\tau -s)\lambda _k(\theta )} u(s) {\text {d}}\!s , \quad k=1,...,n. \end{aligned}$$

Step 2: Take piecewise constant inputs

Let $\tau >0$ (not yet specified) and consider the following partition

$$\begin{aligned}{}[0,N\tau ] = [0,\tau ) \cup [\tau , 2\tau ) \cup \cdots \cup [(N-1)\tau , N\tau ], \end{aligned}$$

where every interval $I_l:=[l\tau , (l+1)\tau )$ has length $\tau$. Further, we consider inputs $u:[0,T] \rightarrow {\mathbb {C}}$ that are constant on every interval, i.e.,

$$\begin{aligned} u\vert _{I_l}(t):= u_l \in {\mathbb {C}}, \quad l =0,...,N-1. \end{aligned}$$

(34)

Let ${\textbf{1}}_{I_l}$ denote the indicator function defined by

$$\begin{aligned} s \mapsto {\textbf{1}}_{I_l}(s) = {\left\{ \begin{array}{ll} 1 &{} \text { if } s \in I_l\\ 0 &{} \text { else }. \end{array}\right. } \end{aligned}$$

Then we take input functions u of the form

$$\begin{aligned} u(T-s)= \sum _{l=0}^{N-1} u_l \, {\textbf{1}}_{I_l}(T-s)= \sum _{l=0}^{N-1} u_{N-1-l} \, {\textbf{1}}_{I_l}(s). \end{aligned}$$

If $\lambda _k(\theta )\ne 0$ the kth component of the solution is then

$$\begin{aligned} \begin{aligned} {\tilde{\varphi }}_k(N\tau ,u,0)(\theta )&= \int _0^{N\tau } {\text {e}}^{ \lambda _k(\theta )s} u(N\tau -s) {\text {d}}\!s \\&= \sum _{l=0}^{N-1} \int _{l\tau }^{(l+1)\tau } \!\!{\text {e}}^{\lambda _k(\theta )s} u_{N-1-l} {\textbf{1}}_{I_l}(s) {\text {d}}\!s\\&= \left( \frac{ {\text {e}}^{\tau \lambda _k(\theta )} -1}{\tau \lambda (\theta )}\right) \sum _{l=0}^{N-1} \tau u_{N-1-l} {\text {e}}^{l \tau \lambda _k(\theta )} . \end{aligned} \end{aligned}$$

If $\lambda _k(\theta )=0$, we have

$$\begin{aligned} {\tilde{\varphi }}_k(N\tau ,u,0)= \sum _{l=0}^{N-1} \int _{l\tau }^{(l+1)\tau } \!\!u_{N-1-l} {\textbf{1}}_{I_l}(l\tau + s) {\text {d}}\!s = \sum _{l=0}^{N-1} \tau \,u_{N-1-l} . \end{aligned}$$

Step 3: Approximate the solution by a polynomial

Observing that

$$\begin{aligned} \sum _{l=0}^{N-1} u_{N-1-l} \left( {\text {e}}^{\tau \lambda _k(\theta )}\right) ^l \end{aligned}$$

defines a polynomial whose coefficients are given by the values of the input function u defined in (34), we define

$$\begin{aligned} p(z) := \sum _{l=0}^{N-1} u_{N-1-l} \, z^l. \end{aligned}$$

(35)

In terms of the polynomial p, the kth component of the solution reads as

$$\begin{aligned} {\tilde{\varphi }}_k(N\tau ,u,0)= {\left\{ \begin{array}{ll} \left( \frac{ {\text {e}}^{\tau \lambda _k(\theta )} -1}{\tau \lambda _k(\theta )}\right) \,\tau \, p\left( {\text {e}}^{\tau \lambda _k(\theta )}\right) &{} \quad \text { if }\lambda _k(\theta )\ne 0\\ \tau \, \, p\left( 1\right) &{} \quad \text { if } \lambda _k(\theta )=0. \end{array}\right. } \end{aligned}$$

(36)

Moreover, since

$$\begin{aligned} \lim _{\tau \rightarrow 0} \frac{{\text {e}}^{\tau z} -1}{ \tau z} = 1 \quad \text { for all } z \in {\mathbb {C}} \setminus \{ 0\}, \end{aligned}$$

there is a $\tau _1>0$ so that for any $\tau \in (0,\tau _1)$ we have

$$\begin{aligned} \left| \frac{{\text {e}}^{\tau \lambda _k(\theta )} -1}{ \tau \lambda _k(\theta )} -1\right| < \tfrac{\varepsilon }{2} \qquad \text { for all }\, \, \lambda _k(\theta )\ne 0, \, \, k=1,...,n. \end{aligned}$$

(37)

Then, by (36) and (37), for every $\tau \in (0, \tau _1)$ each component of the solution satisfies

$$\begin{aligned} \left| \tau p({\text {e}}^{\tau \lambda _k(\theta )})-\varphi _k(N\tau ,u,0)(\theta ) \right| < \tfrac{\varepsilon }{2}. \end{aligned}$$

(38)

The significance of (38) is that it is independent of the number of input values $u_0,...,u_{N-1}$.

Step 4: Approximation of each component of the target states by a polynomial

First, we note that there is a $\tau _2>0$ such that the mappings

$$\begin{aligned} \theta \mapsto {\text {e}}^{\tau \, \lambda _k(\theta )} \end{aligned}$$

are injective for all $k=1,...,n$ and all $\tau \in (0,\tau _2)$. Thus, we fix some $\tau <\min \{\tau _1,\tau _2\}$. The sets

$$\begin{aligned} \varOmega _k:=\{ {\text {e}}^{\tau \lambda _k(\theta )}: \theta \in {\textbf{P}}\} \subset {\mathbb {C}} \end{aligned}$$

define Jordan arcs. Moreover, let $g(\theta ) =T(\theta )f(\theta )$ denote the transformed family of terminal states, we consider the continuous function

$$\begin{aligned} {\tilde{g}}_k :\varOmega _k \rightarrow {\mathbb {C}}, \quad {\tilde{g}}_k(z) = \tfrac{1}{\tau } g_k \left( \lambda ^{-1}_k \left( \tfrac{\ln z}{\tau }\right) \right) . \end{aligned}$$

Applying the results from Sect. 2, we get polynomials $p_1,...,p_n$ such that for each $k=1,...,n$ we have

$$\begin{aligned} |{\tilde{g}}_k(z) - p_k(z)| < \tfrac{\varepsilon }{6\,\tau \, \Vert T\Vert _{M,\infty }} \qquad \forall \, \,z \in \varOmega _k. \end{aligned}$$

Note that, the latter is equivalent to

$$\begin{aligned} | g_k(\theta ) - \tau p_k({\text {e}}^{\tau \lambda _k(\theta )} )| < \tfrac{\varepsilon }{6\, \Vert T\Vert _{M,\infty }} \qquad \forall \, \,\theta \in {\textbf{P}}. \end{aligned}$$

(39)

Step 5: Construction of a single polynomial

Let $\varOmega = \bigcup _{k=1}^n \varOmega _k$. Note that the assumptions imply that the sets $\varOmega _1,..., \varOmega _n$ are disjoint. Consider the holomorphic functions $h_1,...,h_n$ defined by

$$\begin{aligned} h_k:\varOmega \rightarrow {\mathbb {C}}, \quad h_k(z)= {\left\{ \begin{array}{ll} 1 &{} \text { if } z \in \varOmega _k \\ 0 &{} \text { if } z \in \varOmega \setminus \varOmega _k \end{array}\right. } \end{aligned}$$

and compute via Runge’s little theorem polynomials $q_1,...,q_n$ such that

$$\begin{aligned} \sup _{z \in \varOmega } | q_k(z) - h_k(z))| < \frac{\varepsilon }{6\, \tau \, \Vert T\Vert _{M,\infty } \, \sum _{l=1}^n \alpha _{k,l}} \end{aligned}$$

for all $k=1,...,n$, where $\alpha _{k,l}:= \sup _{\theta \in {\textbf{P}}}|p_l(\textrm{e}^{\tau \lambda _k(\theta )})|$. Then we define the polynomial

$$\begin{aligned} p(z) = \sum _{k=1}^n p_k(z)q_k(z). \end{aligned}$$

(40)

Step 6: Definition of the input values

Determine a monomial representation of the polynomial p in (40), i.e.,

$$\begin{aligned} p(z) = \sum _{k=1}^n p_k(z)q_k(z) = \sum _{k=0}^K c_k z^k \end{aligned}$$

(41)

and set $N=K+1$ and define the input $u :[0,N\tau ]\rightarrow {\mathbb {C}}$ by

$$\begin{aligned} u|_{I_l} (t) := c_l, \quad l=0,...,N-1. \end{aligned}$$

Then, putting things together, for given $f \in C_n({\textbf{P}})$ and $\varepsilon >0$ we get for each component

$$\begin{aligned} |\varphi _k(N\tau ,u,0)(\theta )- f_k(\theta )|&\le \Vert T\Vert _{M,\infty } \, \left| {\tilde{\varphi }}_k(N\tau ,u,0)(\theta )-\tau p({\text {e}}^{\tau \lambda _k(\theta )} ) \right| \\ {}&\quad + \Vert T\Vert _{M,\infty } \, \left| \tau p({\text {e}}^{\tau \lambda _k(\theta )} )- g_k(\theta ) \right| \\&< \tfrac{\varepsilon }{2} + \Vert T\Vert _{M,\infty }\, \left| g_k(\theta ) - \tau p_k({\text {e}}^{\tau \lambda _k(\theta )} ) q_k( {\text {e}}^{\tau \lambda _k(\theta )}) \right| \\&\quad + \Vert T\Vert _{M,\infty } \,\tau \,\sum _{l\ne k} \left| p_l({\text {e}}^{\tau \lambda _k(\theta )}) \right| \,\left| q_l({\text {e}}^{\tau \lambda _k(\theta )}) \right| \\&< \tfrac{\varepsilon }{2} + \tfrac{\varepsilon }{6} + \Vert T\Vert _{M,\infty } \, \left| g_k(\theta ) - \tau p_k({\text {e}}^{\tau \lambda _k(\theta )} ) \right| \\&\quad + \Vert T\Vert _{M,\infty } \, \tau \left| 1-q_k( {\text {e}}^{\tau \lambda _k(\theta )})\right| \, \left| p_k({\text {e}}^{\tau \lambda _k(\theta )} ) \right| \\&< \tfrac{\varepsilon }{2} + \tfrac{\varepsilon }{6} + \tfrac{\varepsilon }{6} + \tfrac{\varepsilon }{6} = \varepsilon . \end{aligned}$$

Consequently, we have

$$\begin{aligned} \Vert f - \varphi (N\tau ,u,0)\Vert _{\infty } < \varepsilon . \end{aligned}$$

Remark 2

The approach in Sect. 4 was used in a first and direct proof of the sufficiency of (N1), (N2) and (S2) for uniform ensemble reachability for single-input systems, cf. [28]. Indeed, since the set $\varOmega$ is compact with empty interior and ${\mathbb {C}}\setminus \varOmega$ is connected, it is possible to applying Mergelyan’s theorem directly to the continuous function

$$\begin{aligned} g :\varOmega \rightarrow {\mathbb {C}}, \quad g{|_{\varOmega _k}}(z) = \tfrac{1}{\tau } g_k \left( \lambda ^{-1}_k \left( \tfrac{\ln z}{\tau }\right) \right) . \end{aligned}$$

Thus, there is a polynomial p such that

$$\begin{aligned} |g(z) - \tau p(z)|< \varepsilon \qquad \forall \, z \in \varOmega . \end{aligned}$$

Also we note that the property that the interior of $\varOmega$ is empty is a special case of Mergelyan’s result, which was proven earlier by Lavrientev in 1934, cf. [24]. However, Lavrientev’s proof is not available to the author. Maybe the proof can be used to obtain another constructive method, at least for special cases.

5 Comments

In this section we discuss some perspectives on the methods presented in Sect. 3.

Degrees of freedom

If the approximation domains are given by Jordan arcs, the methods for (S1) and (S2) have some degrees of freedom. For instance, there are plenty of ways to close a Jordan arc to a Jordan curve. This in turn frees up space for designing the conformal mappings such that conditions (21), (23) and (28), (30) are satisfied. Thus, the computation of the conformal mappings and the design of the Jordan curves can be treated simultaneously. Besides, for pairs (A, b) where A is block diagonal and each block satisfies the sufficiency condition (S1) or (S2), the methods can be applied for every block and an overall polynomial can be constructed in the same manner the polynomials $q_1,...,q_n$ are determined in the methods for (S2).

Real approximation domains

In case the domains, where the approximations take place, are compact real intervals, the construction process is much easier. This can be achieved by considering a mixture of open-loop and feedback control inputs of the form

$$\begin{aligned} K(\theta )x_t(\theta ) + u_t. \end{aligned}$$

Indeed, due to the necessary condition (N1), a suitable feedback term $K(\theta )x_t(\theta )$ can be used to shift the spectra of $A(\theta )$ to an appropriate position. That is, via $K(\theta )$ the matrices $A(\theta )+B(\theta )K(\theta )$ can have any spectra. Equivalently, the coefficients of the characteristic polynomials of $A(\theta )+B(\theta )K(\theta )$ can be arbitrarily designed. For more details on this we refer to [20, 44, Section 6.4].

Multi-input systems

The known necessary and sufficient conditions for ensemble reachability are in the multi-input case much less precise compared to the single-input systems. For special classes of systems, such as $(A(\theta ),B(\theta ))$ upper triangular, the problem can be traced back to the single-input case. In this context, a general but very strong sufficiency condition is that the Hermite indices of the pair $(A(\theta ),B(\theta ))$ are constant. Technicalities aside, these systems can be transformed into the Hermite canonical form

$$\begin{aligned} \begin{pmatrix} {\tilde{A}}_{11}(\theta ) &{}\quad \cdots &{}\quad {\tilde{A}}_{1k}(\theta ) \\ &{}\quad \ddots &{}\quad \vdots \\ 0 &{}\quad &{}\quad {\tilde{A}}_{kk}(\theta ) \end{pmatrix}, \quad \begin{pmatrix} {\tilde{b}}_{1}(\theta ) &{}\quad &{}\quad 0 &{}\quad * &{}\quad * \\ &{}\quad \ddots &{}\quad &{}\quad *&{}\quad * \\ 0 &{}\quad &{}\quad {\tilde{b}}_{k}(\theta )&{}\quad *&{}\quad * \end{pmatrix}, \end{aligned}$$

and the construction problem can be solved by applying the single-input methods to the subpairs $({\tilde{A}}_{ii}(\theta ),{\tilde{b}}_i(\theta ))$. We note that it is not required to use one of the methods for all subpairs. In particular, this applies to the often studied class $(\theta A,B)$. For more details, we refer to [14] and the references therein.