An Iterative Wiener--Hopf method for triangular matrix functions with exponential factors

TL;DR

This paper presents an iterative Wiener--Hopf method for approximating solutions to a class of equations involving triangular matrix functions with exponential factors, enabling controlled error reduction and practical application in acoustics.

Contribution

It introduces a novel iterative approach generalizing pole removal techniques for Wiener--Hopf equations with explicit error estimates and efficient convergence.

Findings

Error can be explicitly estimated and minimized.

Few iterations suffice for practical accuracy.

Method successfully applied to acoustics problems.

Abstract

This paper introduces a new method for constructing approximate solutions to a class of Wiener--Hopf equations. This is particularly useful since exact solutions of this class of Wiener--Hopf equations, at the moment, cannot be obtained. The proposed method could be considered as a generalisation of the pole removal technique. The error in the approximation can be explicitly estimated, and by a sufficient number of iterations could be made arbitrary small. Typically only a few iterations are required for practical purposes. The theory is illustrated by numerical examples that demonstrate the advantages of the proposed procedure. This method was motivated and successfully applied to problems in acoustics.

Equations

$$\left(\begin{array}{c}{\Phi }_{-}^{(0)}(\alpha )\\ {\Phi }_{-}^{(L)}(\alpha )\end{array}\right)=\left(\begin{array}{cc}A(\alpha )& B(\alpha )e^{i\alpha L}\\ C(\alpha )e^{-i\alpha L}& 0\end{array}\right)\left(\begin{array}{c}{\Psi }_{+}^{(0)}(\alpha )\\ {\Psi }_{+}^{(L)}(\alpha )\end{array}\right)+\left(\begin{array}{c}{f}_1(\alpha )\\ f_2(\alpha )\end{array}\right),$$

$$|F(x_1)-F(x_2)|\le C|x_1-x_2{|}^{\lambda },\text{ for all }{x}_1,x_2\in \mathbb{R}$$

$$|F(x_1)-F(x_2)|\le C{\left|\frac{1}{x_1}-\frac{1}{x_2}\right|}^{\lambda },$$

$$\{\{[0]\}\}=L_2(\mathbb{R})\cap \text{ Ho¨lder},$$

$$f_{+}(t)=\left\{\begin{array}[]{rl}f(t)&\text{if }t>0,\\ 0&\text{if }t<0,\\ \end{array}\right.\quad f_{-}(t)=\left\{\begin{array}[]{rl}0&\text{if }t>0,\\ -f(t)&\text{if }t<0.\\ \end{array}\right.$$

$$F_{+}(z)\in \{\{[a,\infty ]\}\}$$

$$F_{-}(z)\in \{\{[-\infty ,b]\}\}$$

$$F(z)\in \{\{[a,b]\}\}$$

$$F(t)=F^{+}(t)+F^{-}(t),$$

$$K(t)=K_{+}(t)K_{-}(t),$$

$$H(z)\in \{\{[\text{max}(a,\alpha ),\text{min}(b,\beta )]\}\},$$

$$K(t){\Phi }_{-}(t)={\Phi }_{+}(t)+F(t).$$

$$\Phi (\alpha ,y)=$$

$$=$$

$${\Phi }^{(L)}(\alpha ,y)=$$

$$=$$

$${\Phi }^{(L)}(\alpha ,y)=\Phi (\alpha ,y)e^{i\alpha L}.$$

$$\left(\begin{array}{c}{\Phi }_{-}^{(0)}(\alpha )\\ {\Phi }_{-}^{(L)}(\alpha )\end{array}\right)=-\left(\begin{array}{cc}\frac{1-\gamma (\alpha )P}{\gamma (\alpha )}& P{e}^{i\alpha L}\\ \frac{1}{\gamma (\alpha )}{e}^{-i\alpha L}& 0\end{array}\right)\left(\begin{array}{c}{\Phi }_{+}^{{}^{\prime }(0)}(\alpha )\\ {\Phi }_{+}^{{}^{\prime }(L)}(\alpha )\end{array}\right)+\left(\begin{array}{c}{f}_1(\alpha )\\ f_2(\alpha )\end{array}\right),$$

$$\left(\begin{array}{c}{\Phi }_{-}^{{}^{\prime }(0)}(\alpha )\\ {\Phi }_{-}^{(L)}(\alpha )\end{array}\right)=-\left(\begin{array}{cc}\gamma (\alpha )& e^{i\alpha L}\\ -e^{-i\alpha L}& 0\end{array}\right)\left(\begin{array}{c}{\Phi }_{+}^{(0)}(\alpha )\\ {\Phi }_{+}^{{}^{\prime }(L)}(\alpha )\end{array}\right)-\left(\begin{array}{c}{g}_1(\alpha )\\ 0\end{array}\right),$$

$$\left(\begin{array}{cc}\frac{-e^{-i\alpha L}}{B_{-}(\alpha )}& \frac{A(\alpha )}{C(\alpha )B_{-}(\alpha )}\\ \frac{1}{B_{-}(\alpha )}& 0\end{array}\right)\left(\begin{array}{c}{\Phi }_{-}^{(0)}(\alpha )\\ {\Phi }_{-}^{(L)}(\alpha )\end{array}\right)$$

$$=\left(\begin{array}{cc}0& -B_{+}(\alpha )\\ \frac{A(\alpha )}{B_{-}(\alpha )}& B_{+}(\alpha )e^{i\alpha L}\end{array}\right)\left(\begin{array}{c}{\Psi }_{+}^{(0)}(\alpha )\\ {\Psi }_{+}^{(L)}(\alpha )\end{array}\right)+\left(\begin{array}{c}{f}_3\\ f_4\end{array}\right),$$

$$f_3=\frac{-e^{-i\alpha L}}{B_{-}(\alpha )}{f}_1+\frac{A(\alpha )}{C(\alpha )B_{-}(\alpha )}{f}_2\text{and}{f}_4=\frac{f_1}{B_{-}(\alpha )}.$$

$$\left(\begin{array}{cc}k{B}_{-}(\alpha )C_{+}(\alpha )& B(\alpha )e^{i\alpha L}\\ C(\alpha )e^{-i\alpha L}& 0\end{array}\right),$$

$$\left(\begin{array}{cc}\frac{-e^{-i\alpha L}}{B_{-}(\alpha )}& {\left(\frac{A(\alpha )}{C(\alpha )B_{-}(\alpha )}\right)}^{-}\\ \frac{1}{B_{-}(\alpha )}& 0\end{array}\right)\left(\begin{array}{c}{\Phi }_{-}^{(0)}(\alpha )\\ {\Phi }_{-}^{(L)}(\alpha )\end{array}\right)+\left(\begin{array}{c}{\left(\frac{A(\alpha )}{C(\alpha )B_{-}(\alpha )}\right)}^{+}{\Phi }_{-}^{(L)}(\alpha )\\ -{\left(\frac{A(\alpha )}{B_{-}(\alpha )}\right)}^{-}{\Psi }_{+}^{(0)}(\alpha )\end{array}\right)$$

$$=\left(\begin{array}{cc}0& -B_{+}(\alpha )\\ {\left(\frac{A(\alpha )}{B_{-}(\alpha )}\right)}^{+}& B_{+}(\alpha )e^{i\alpha L}\end{array}\right)\left(\begin{array}{c}{\Psi }_{+}^{(0)}(\alpha )\\ {\Psi }_{+}^{(L)}(\alpha )\end{array}\right)+\left(\begin{array}{c}{f}_3\\ f_4\end{array}\right),$$

$$\frac{-e^{-i\alpha L}}{B_{-}}{\Phi }_{-}^{(0)}+{\left(\frac{A}{C{B}_{-}}\right)}^{-}{\Phi }_{-}^{(L)}+{\left({\left(\frac{A}{C{B}_{-}}\right)}^{+}{\Phi }_{-}^{(L)}\right)}^{-}-f_3^{-}=$$

$$-B_{+}{\Psi }_{+}^{(L)}-{\left({\left(\frac{A}{C{B}_{-}}\right)}^{+}{\Phi }_{-}^{(L)}\right)}^{+}+f_3^{+}=$$

$$\frac{1}{B_{-}}{\Phi }_{-}^{(0)}-{\left({\left(\frac{A}{B_{-}}\right)}^{-}{\Psi }_{+}^{(0)}\right)}^{-}-f_4^{-}=$$

$${\left(\frac{A}{B_{-}}\right)}^{+}{\Psi }_{+}^{(0)}+B_{+}{e}^{i\alpha L}{\Psi }_{+}^{(L)}+{\left({\left(\frac{A}{B_{-}}\right)}^{-}{\Psi }_{+}^{(0)}\right)}^{-}+f_4^{+}=$$

$${\left(\frac{A}{C{B}_{-}}\right)}^{-}{\Phi }_{-}^{(L)}=$$

$$B_{+}{\Psi }_{+}^{(L)}=$$

$${\left(\frac{A}{B_{-}}\right)}^{+}{\Psi }_{+}^{(0)}=$$

$$\frac{{\Phi }_{-}^{(0)}}{B_{-}}=$$

$$K_1^{-}{\Phi }_{-}^{(L)}=$$

$$K_2^{+}{\Psi }_{+}^{(L)}=$$

$$K_3^{+}{\Psi }_{+}^{(0)}=$$

$$K_4^{-}{\Phi }_{-}^{(0)}=$$

$$K_1^{-}{\Phi }_{-}^{(L)}=$$

Full text

1. Introduction

Many boundary value problems in mathematical physics can be approached by the Wiener–Hopf method. Originally the Wiener–Hopf technique was developed for linear PDEs with semi-infinite boundary conditions like the Sommerfeld half-plane problem. This was later extended to boundary conditions on multiple semi-infinite lines [history]. Although, the reduction to the Wiener–Hopf equation is still straightforward, finding a solution of the resulting equation became challenging [Constructive_review, Ant_crack, Mer, Daniele, Fokas_uni, Crowdy-Elena]. This is due to the fact that the Wiener–Hopf factorisation of a matrix (rather than a scalar) function is now needed. Hence the ability to solve such equations is crucial to extending the classical use of Wiener–Hopf techniques to more realistic and complicated settings. Also, such types of matrix Wiener–Hopf equations are associated with convolution-type operators on a finite interval [Feldman2000, Feldman2004, Pure_rew_2013] and arise in a number of applications [Nigel_cyl, Abr_exp]. The aim of this paper is to develop an algorithmic iterative method of solution for some equations of this type.

More precisely, we construct an approximate solution of a Wiener–Hopf equation with triangular matrix functions containing exponential factors. The aim is to find functions $\Phi_{-}^{(0)}(\alpha)$, $\Phi_{-}^{(L)}(\alpha)$, $\Psi_{+}^{(0)}(\alpha)$ and $\Psi_{+}^{(L)}(\alpha)$ analytic in respective half-planes, satisfying the following relationship

[TABLE]

on the strip $a\leq\Im(\alpha)\leq b$. The remaining functions $A(\alpha)$, $B(\alpha)$ and $C(\alpha)$ are known and $L$ is a positive constant. The conditions on the matrix functions are specified in Section 2. The existence of a such factorisation under certain assumptions was addressed in [Spit_book, p. 150] and [Existance_osi].

Wiener–Hopf equations of the type (1) have been the topic of previous research, for example, in the case of meromorphic matrix entries [Antipov, Aktosun1992] and in the framework of almost periodic functions [karlovich-spit-alm-per]. However, presently no complete solution of (1) is known. Furthermore, the general question of constructive Wiener–Hopf factorisation is widely open [history, Constructive_review].

In the literature, concerning applications of the Wiener–Hopf technique, one of most widely used methods is the so-called “pole removal” or “singularities matching” [bookWH]§ 4.4, 5.3 [daniele2014wiener]§ 4.4.2. It has a severe limitation that certain functions have to be rational or meromorphic. One way to extend the use of this method is by employing a rational approximation [Pade, My1], which was successfully used in [Abraha_all_pade, Ab_ex, Nigel_cyl, MyD-K]. However, even with this extension the class of functions, which can be solved, is rather limited. In this paper we propose a different extension to the pole removal technique for functions that have arbitrary singularities. This work was motivated by certain problems in acoustics which are discussed below.

We note, that there is an additional difficulty in finding a factorisation of the matrix in (1) because of the presence of analytic functions $e^{i\alpha L}$ and $e^{-i\alpha L}$ which have exponential growth in one of the half-planes. The first step of the procedure proposed here is a partial factorisation, that has at most polynomial growth in the respective half-planes. The classes for which this is a complete factorisation are also discussed. The next step is the Wiener–Hopf additive splitting of the remainder term that hinders the application of Liouville’s theorem. The additive Wiener–Hopf splitting is routine unlike the multiplicative one, which is also utilised in other novel methods [Mishuris-Rogosin, assym_2]. After the application of the analytic continuation there are still some unknowns in the formula, those are approximated by an iterative procedure. The presence of the exponential terms speeds up convergence. At each step a scalar Wiener–Hopf equation is solved. We will compare the proposed method to the pole removal technique. It is shown that the iterations converge quickly to the exact solution.

The procedure could be summarised as follows (see Section 3 for details):

1. (1)

A partial factorisation with exponential factors in the desired half-planes. 2. (2)

Additive splitting of some terms. 3. (3)

Application of Liouville’s theorem. 4. (4)

Iterative procedure to determine the remaining unknowns.

This procedure bypasses the need to construct a multiplicative matrix factorisation. So in particular partial indices (which are know to be linked to stability [Spit_book, MyD-K]) are not obtained. Instead the growth at infinity of certain terms play a role. In this paper we will treat the base case with no growth at infinity and some other cases will be treated in [my+Lorna].

The structure of the paper is as follows. In Section 2 the required classes of functions are introduced and their essential properties are listed. We also provide some motivation behind those Wiener–Hopf systems. In Section 3 the proposed iterative procedure is described in detail and its convergence is examined in Section 4. Section 5 presents numerical results of two examples (graphically illustrated) to compare the iterative procedure to the exact solution in a variety of cases. Lastly we describe possible future work.

2. Preliminaries

In order to formulate the problem (1) we have to specify the suitable class of functions for the all terms in the equation. Those classes of functions are convenient for formulating scalar Wiener–Hopf splittings which will form the foundation of the proposed method. We will also need some properties of this class of function which will ensure we stay the the desired class of functions after each iterations. We will also briefly describe the motivation behind this method. Firstly, we will recall some definitions.

2.1. Classes of functions $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$

The Hölder continuous functions on the compactified real line are defined as functions $F(x)$ for which there exist constants $C$ and $\lambda$ such that for all real $x_{1}$ and $x_{2}$ we have

[TABLE]

and for all real $x_{1}$ and $x_{2}$ with modulus greater than one the following holds

[TABLE]

which is a Hölder condition around infinity. In the following, the “real line” will always mean the “compactified real line”.

A Hölder continuous function produces well-defined boundary values $G_{\pm}(t)=\lim_{y\searrow 0}G(t\pm iy)$ of its Cauchy type integral $G(t+iy)$ at every point $t$ of the compactified real line [Musk]. The class $L_{2}(\mathbb{R)}$ is also very useful since the Fourier transform is an isometry of $L_{2}(\mathbb{R)}$ due to Plancherel’s theorem. Thus, the intersection [Ga-Che]*§ 1.2:

[TABLE]

turns out to be convenient for the Wiener–Hopf problems. The pre-image of $\{\kern-3.00003pt\{\kern-2.5pt[0]\kern-2.5pt\}\kern-3.00003pt\}$ under the Fourier transform is denoted $\{\kern-2.5pt[0]\kern-2.5pt\}$.

Given a function $f\in\{\kern-2.5pt[0]\kern-2.5pt\}$ on the real line we can define the splitting

[TABLE]

Using this splitting, we define the class $\{\kern-2.5pt[0,\infty]\kern-2.5pt\}$ to contain functions $f_{+}(t)$ and the class $\{\kern-2.5pt[-\infty,0]\kern-2.5pt\}$ to contain functions $f_{-}(t)$ for all $f\in\{\kern-2.5pt[0]\kern-2.5pt\}$.

In the rest of the paper we will need to refer to functions that are analytic on strips or (shifted) half-planes. Following [Ga-Che]*§ 13, we define $f\in\{\kern-2.5pt[a]\kern-2.5pt\}$ if $e^{-ax}f\in\{\kern-2.5pt[0]\kern-2.5pt\}$, that is a shift in the Fourier space. Finally, $f=f_{+}+f_{-}\in\{\kern-2.5pt[a,b]\kern-2.5pt\}$ if $f_{+}\in\{\kern-2.5pt[a]\kern-2.5pt\}$ and $f_{-}\in\{\kern-2.5pt[b]\kern-2.5pt\}$. From the definition of $f_{+}$ and $f_{-}$ it is clear that also $f_{+}\in\{\kern-2.5pt[a,\infty]\kern-2.5pt\}$ and $f_{-}\in\{\kern-2.5pt[-\infty,b]\kern-2.5pt\}$.

The Fourier transform of functions in the class $\{\kern-2.5pt[a,b]\kern-2.5pt\}$ is denoted $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$ (including the case $a=-\infty$ and/or $b=\infty$). The celebrated Paley–Wiener theorem states that the following

[TABLE]

For the remainder of the paper we will assume that $a<0$ and $b>0$, i.e. the strip encloses the real axis.

Now we will recall the additive and multiplicative scalar splitting of functions belonging to $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$. A convention will be used to distinguish the additive and the multiplicative splitting by using the superscript and subscript notations: e.g $F^{\pm}$ for additive and $K_{\pm}$ for multiplicative ones.

Theorem 2.1** (Additive splitting [Ga-Che]).**

On the real line a function $F(t)\in\{\kern-2.5pt[0]\kern-2.5pt\}$ is given. There exist two functions $F^{\pm}(z)$ analytic in the upper and lower half-planes with boundary functions on the real line belonging to the classes $\{\kern-3.00003pt\{\kern-2.5pt[0,\infty]\kern-2.5pt\}\kern-3.00003pt\}$ and $\{\kern-3.00003pt\{\kern-2.5pt[-\infty,0]\kern-2.5pt\}\kern-3.00003pt\}$, satisfying:

[TABLE]

on the real line.

Next, the multiplicative splitting or factorisation problem is examined. The index of a continuous non-zero function $K(t)$ on the real line is the winding number of the curve $(\text{Re }K(t),\mathrm{Im}\,K(t))$, $t\in\mathbb{R}$.

Theorem 2.2** (Multiplicative splitting [Ga-Che]).**

Let a non-zero function $K(t)$, such that $K(t)-1\in\{\kern-2.5pt[0]\kern-2.5pt\}$ and $\mathrm{ind}\,K(t)=0$, be given. There exist two functions $K_{\pm}(z)$ analytic in the upper and lower half-planes respectively with boundary functions $K_{\pm}(t)-1$ on the real line belonging to the classes $\{\kern-3.00003pt\{\kern-2.5pt[0,\infty]\kern-2.5pt\}\kern-3.00003pt\}$ and $\{\kern-3.00003pt\{\kern-2.5pt[-\infty,0]\kern-2.5pt\}\kern-3.00003pt\}$, satisfying:

[TABLE]

on the real line.

Remark 2.3*.*

Note, that if the original function in Theorem 2.1 or Theorem 2.2 is from the class $\{\kern-2.5pt[a]\kern-2.5pt\}$ then the components belong to the classes $\{\kern-3.00003pt\{\kern-2.5pt[a,\infty]\kern-2.5pt\}\kern-3.00003pt\}$ and $\{\kern-3.00003pt\{\kern-2.5pt[-\infty,a]\kern-2.5pt\}\kern-3.00003pt\}$.

With those definitions we can specify the classes of functions in equation (1). The given functions $A(\alpha)-1$, $C(\alpha)-1$, $B(\alpha)-1$ are required to be in $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$ and $A(\alpha)$, $C(\alpha)$, $B(\alpha)$ to have zero index. We also need that $C(\alpha)$ and $B(\alpha)$ have no zeroes on the strip—this corresponds to the determinant of the original Wiener–Hopf matrix (1) being non-zero. We will also require that $A(\alpha)$ has no zeroes on the strip, this corresponds to the sub-problem of (1) (when $L$ is very large) being non-degenerate. The “forcing terms” $f_{1}(\alpha)$ and $f_{2}(\alpha)$ should be in $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$. We look for $\Phi_{-}^{(0)}(\alpha)$ and $\Phi_{-}^{(L)}(\alpha)$ in $\{\kern-3.00003pt\{\kern-2.5pt[-\infty,b]\kern-2.5pt\}\kern-3.00003pt\}$ and $\Psi_{+}^{(0)}(\alpha)$ and $\Psi_{+}^{(L)}(\alpha)$ in $\{\kern-3.00003pt\{\kern-2.5pt[a,\infty]\kern-2.5pt\}\kern-3.00003pt\}$.

2.2. Properties of functions in $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$

In this subsection we briefly describe some of the properties of $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$ which will ensure that the solution we seek is in the desired class. We will need the fact that if $G\in\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$, $F\in\{\kern-3.00003pt\{\kern-2.5pt[\alpha,\beta]\kern-2.5pt\}\kern-3.00003pt\}$ ($a<b$, $\alpha<\beta$) and $H(z)=F(z)G(z)$ then

[TABLE]

i.e $H(z)$ is analytic in the respective strip. This is a simple but still fundamental property for the rest of the derivation. Note, that a similar statement does not hold for the class of $L_{2}(\mathbb{R})$.

We will also use that if $G-1\in\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$, $\mathrm{ind}\,G(t)=0$ and non-zero on the strip $a\leq\Im(\alpha)\leq b$ then $G^{-1}-1\in\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$.

Note that in the proposed method we only solve scalar Wiener-Hopf equations of the similar form to the ones detailed in [Ga-Che]*§ 3.2. That is for a scalar equations on the real line

[TABLE]

we need that $K$ is non-zero function, such that $K(t)-1\in\{\kern-2.5pt[0]\kern-2.5pt\}$ and $\mathrm{ind}\,K(t)=0$. Also require that $F(t)\in\{\kern-2.5pt[0]\kern-2.5pt\}$ then we will can find $\Phi_{+}$ and $\Phi_{-}$ in $(\{\kern-3.00003pt\{\kern-2.5pt[0,\infty]\kern-2.5pt\}\kern-3.00003pt\}$ and $\{\kern-3.00003pt\{\kern-2.5pt[-\infty,0]\kern-2.5pt\}\kern-3.00003pt\}$ on the real axis. We will need the extension of those result for the case when (6) holds on a strip $a\leq\Im(\alpha)\leq b$. This has been obtained in [MyStrip].

2.3. Motivation

The initial motivation for this work came from the following acoustics problem: investigate the effect of a finite poroelastic trailing edge on noise production [owl, Lorna]. The situation is modelled with a plane-wave scattering by a rigid half-line plate $(-\infty,0)$ and an poroelastic edge $(0,L)$ with the poroelastic to rigid plate transition at $x=0$. The matrix Wiener–Hopf problem is obtained from consideration of the Fourier transform with respect to the ends of the poroelastic plate at $x=0$ and $x=L$. Define the half-range and full-range Fourier transform with respect to $x=0$:

[TABLE]

In the case when there is only one transition point $x=0$, the unknown functions would be $\Phi^{(0)}_{-}(\alpha,y)$, $\Phi^{(0)}_{+}(\alpha,y)$ (or their derivatives) and we would only need to solve a scalar Wiener–Hopf equation [bookWH]. Since there is a change of boundary conditions at $x=L$ as well, we will also define the Fourier transforms with respect to the point $x=L$:

[TABLE]

The relation between the transforms is:

[TABLE]

The next step is to write down the relationship between different half-range transforms by using the boundary conditions [daniele2014wiener]. These relations can be combined to form a matrix Wiener–Hopf equation. The resulting matrix, which motivated the present method, is (see [my+Lorna] for a detailed discussion):

[TABLE]

where $\gamma(\alpha)=\sqrt{\alpha^{2}-k^{2}_{0}}$, $k_{0}$ is the acoustic wave number and, in the simplest case, $P$ is a constant. The exponentials are due to the boundary conditions since the half-range Fourier transform with respect to different points are needed.

We note that on the left hand side the unknown functions are minus half-transforms $\Phi_{-}^{(0)}(\alpha)$ and $\Phi_{-}^{(L)}(\alpha)$, and on the right are the derivatives of the plus half-transforms which are in this paper denoted by $\Psi_{+}^{(0)}(\alpha)$ and $\Psi_{+}^{(L)}(\alpha)$. This problem is not treated in this paper further since the purpose of this paper is to present the iterative procedure in the simplest and general case in order to make it easier to apply to many situations. For a detailed discussion of this particular problem see [my+Lorna]. The the later paper we also discuss other examples, the diffraction with finite rigid plate, leading to the following Wiener–Hopf equation.

[TABLE]

The advantage in considering the following example is that although the exact matrix factorisation cannot be constructed, the results of the new iterative procedure can be compared to the exact solution obtained by other methods (Mathieu functions).

3. Iterative Wiener–Hopf Factorisation

The most characteristic aspect of the Wiener–Hopf method is the application of Liouville’s theorem in order to obtain two separate equations from one equation. In order for Liouville’s theorem to be used two conditions have to be satisfied: the analyticity and (at most) polynomial growth at infinity. These two conditions will be treated here in turn. First, a partial factorisation is considered that has the exponential functions in the right place and some of the required analyticity, that is:

[TABLE]

on a strip $a\leq\Im(\alpha)\leq b$, where

[TABLE]

Note, that $f_{3}$ and $f_{4}$ are still in $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$. There are two cases which would allow to solve the above equation exactly. The first case is when $\frac{A(\alpha)}{C(\alpha)B_{-}(\alpha)}$ is in $\{\kern-3.00003pt\{\kern-2.5pt[-\infty,b]\kern-2.5pt\}\kern-3.00003pt\}$ and $\frac{A(\alpha)}{B_{-}(\alpha)}$ is in $\{\kern-3.00003pt\{\kern-2.5pt[a,\infty]\kern-2.5pt\}\kern-3.00003pt\}$. Then the matrix factorisation has been already achieved in 11. In particular, this is true for matrices that have the form:

[TABLE]

where $k$ is a constant. The second important case is when $\left(\frac{A(\alpha)}{C(\alpha)B_{-}(\alpha)}\right)^{+}$ and $\left(\frac{A(\alpha)}{B_{-}(\alpha)}\right)^{-}$ are rational functions. Then the pole removal method [bookWH]§ 4.4, 5.3 [daniele2014wiener]§ 4.4.2 could be employed to obtain the factorisation.

In the generic case we start with the partial factorisation (11) and then use the additive Wiener–Hopf splitting. We present the detailed description now. The equation (11) can be rearranged as

[TABLE]

on a strip $a\leq\Im(\alpha)\leq b$. Recall that the additive splitting is denoted $F^{\pm}$ and the multiplicative splitting by $K_{\pm}$. As the next step we make the additive splittings of the second term of (12), which are possible since it is in the class $\{\kern-3.00003pt\{\kern-2.5pt[a,b]\kern-2.5pt\}\kern-3.00003pt\}$—in the same way as for the second term of (13). Then, the Liouville’s theorem could be applied because the exponential functions are in the correct place and all the functions have desired the analyticity. Thus, we can apply the Wiener–Hopf procedure as usual and the four equations (defined for all $\alpha$ in the complex plane) then become

[TABLE]

Note that Liouville’s theorem is applied before any approximations are made. The four equations can be rearranged:

[TABLE]

When the equations are written in this form it is clear that if $\Phi_{-}^{(L)}$ is known then it could be used to calculate $\Psi_{+}^{(L)}$ and this, in turn, produces $\Psi_{+}^{(0)}$ followed by the calculation of $\Phi_{-}^{(0)}$ and then it loops round. To avoid cumbersome notations, we will label coefficients in (14–17) by $K^{\pm}_{i}$ and obtain the following system:

[TABLE]

Using equations (19) and  (21), we eliminate $\Psi_{+}^{(L)}$ and $\Phi_{-}^{(0)}$ from (18) and (20) respectively:

[TABLE]

Functions $\Psi_{+}^{(L)}$ and $\Phi_{-}^{(0)}$ can be found from (19) and  (21) once (22) and (23) will be solved, where $K_{1}=\frac{A}{CB_{-}}$ and $K_{3}=\frac{A}{B_{-}}$ in the original notation.

So far the equations are exact, but in order to make progress an approximation will be used now. In order to solve approximately we will describe an iterative procedure, where the $n$-th iteration is denoted by $\Phi_{-}^{(L)n}$ and $\Psi_{+}^{(0)n}$. If $\Phi_{-}^{(L)n}$ is known it could be substituted into (23) to calculate $\Psi_{+}^{(0)n+1}$ and then the function $\Psi_{+}^{(0)n+1}$ can be used in (22) to find $\Phi_{-}^{(L)n+1}$ and so on.

Hence, for this iterative procedure it is enough to choose an initial value of $\Phi_{-}^{(L)0}$. Since equation (22) can be considered on a horizontal line $\Im(\alpha)=a<0$, on that line the term with $e^{-i\alpha L}$ will be small (especially for $L$ large). This justifies neglecting the term with $e^{-i\alpha L}$ as a first approximation. Hence to find $\Phi_{-}^{(L)0}$ the aim is to solve

[TABLE]

The above equation can be rearranged as a scalar Wiener–Hopf equation in the following manner

[TABLE]

Introduce an unknown function $D^{+}$ defined by

[TABLE]

Then, combining (25) and (26) we obtain a Wiener–Hopf equation:

[TABLE]

This equation has the form discussed in Section 2.2 so we will have $\Phi_{-}^{(L)0}\in~{}\{\kern-3.00003pt\{\kern-2.5pt[-\infty,a]\kern-2.5pt\}\kern-3.00003pt\}$ as desired. The solution of this equation is

[TABLE]

This can be taken as the initial approximation of the solution, which is used in our iterative procedure.

To compute the next step we will need to solve (23) for $\Psi_{+}^{(0)1}$ on a horizontal line $\Im(\alpha)=b>0$. Note that $\Phi_{-}^{(L)1}$ is defined on this line and can be evaluated numerically. Exactly in the same manner as (24), the solution of

[TABLE]

will lead to the scalar Wiener–Hopf equation now on the line $\Im(\alpha)=b$. It can be found that

[TABLE]

It is easy to see that this formula holds for all iterations with a trivial change of $\Psi_{+}^{(0)1}$ to $\Psi_{+}^{(0)n}$ and $\Phi_{-}^{(L)0}$ to $\Phi_{-}^{(L)n-1}$. Similarly, the general recurrence formula for $\Phi_{-}^{(L)n}$ is

[TABLE]

The convergence of this procedure is examined in the next section. Numerical examples of this procedure are given in Section 5 and are compared to exact solutions (which are known for these special cases).

4. Convergence of the method

The convergence of iterations relies on consideration of equations (22) and (23) on different lines $\Im(\alpha)=a<0$ and $\Im(\alpha)=b>0$ within the strip of analyticity $a\leq\Im(a)\leq b$, in a similar way as it was done in [MyStrip]. We will employ the following notation

[TABLE]

From (28)–(29), the difference of the values of the function after $n+1$ and $n$ times is

[TABLE]

and

[TABLE]

Note that the forcing terms $f_{i}^{\pm}$ do not influence the convergence. In order to estimate the differences (30) and (31) in magnitude we need some inequalities for the Wiener–Hopf additive decomposition. This has been addressed in [My1] for all $L_{p}$ spaces, here we will need a very special case of that result which has a particularly simple form. We will use that if $F(t)=F^{+}(t)+F^{-}(t),$ as in Theorem 2.1, then

[TABLE]

The same inequality can be observed from the fact that the map $F\to F_{\pm}$ in $L_{2}(\mathbb{R})$ is the Szegö orthoprojector with the norm $1$. Since we will be looking at the $L_{2}$ norm on different lines within the strip of analyticity, the $L_{2}$ norm of a function $f(x+ia)$ on the line $\Im(\alpha)=a$ will be denoted by

[TABLE]

As well as employing (30) to assert convergence we will need to derive a relationship for $E_{n+1}^{-}$ and $E_{n}^{-}$. This is done by a similar procedure used to derive (30) but considering the other unknown in the scalar Wiener-Hopf equation. In other words the expression of the unknown plus and minus function are obtained. This is done explicitly below. Write (22) and (23) as they are used on $n$-th iteration:

[TABLE]

We can add $K_{1}^{+}\Phi_{-}^{(L)n}$ to both sides of (33) and $K_{3}^{-}\Psi_{+}^{(0)n}$ to both sides of (34) to obtain

[TABLE]

And if we consider the difference between the consecutive iterations we derive

[TABLE]

This makes the coupling between the equations explicit. Note that

[TABLE]

where $D_{n}^{+}$ can be derived by solving the scalar Wiener–Hopf equation (37) with the unknown minus function $\Phi_{-}^{(L)n}-\Phi_{-}^{(L)n-1}$, unknown plus function $D_{n}^{+}$ and $e^{-i\alpha L}E_{n}^{-}$ is assumed to be known from the previous iterations. Similarly $E_{n+1}^{-}$ is obtained by solving the scalar Wiener–Hopf equation (38).

We will obtain estimates of the size of $\lVert D_{n}^{+}\rVert_{2}^{a}$. Let $\max_{x\in\mathbb{R}}|(K_{1})_{+}(x+ai)|=d_{1}$ and $\max_{x\in\mathbb{R}}|(K_{1})_{+}^{-1}(x+ai)e^{-iL(x+ai)}|=\epsilon_{1}$. Note that $K_{1}$ is bounded since it is in $L_{2}(\mathbb{R})\cap\text{ H\"{o}lder}$ and note that $K_{1}$ is non-zero. Then using (32) we have

[TABLE]

Similarly, defining $\max_{x\in\mathbb{R}}(K_{3})_{-}(x+bi)=d_{2}$ and $\max_{x\in\mathbb{R}}e^{iL(x+bi)}(K_{3})_{-}^{-1}(x+bi)=\epsilon_{2}$, we obtain

[TABLE]

Next, we note that the following is true

[TABLE]

This is intuitively clear since $E_{n+1}^{-}$ is further from singularities on the line $\Im(\alpha)=a$ than on $\Im(\alpha)=b$ and the other way round for $D_{n}^{+}$. The inequalities follows from the Poisson formula for the real line [Nikolski]*Cor 6.4.1 and Hölder inequality. Combining (41) with inequalities (39) and (40) we obtain the key result for demonstrating convergence:

[TABLE]

The convergence of the procedure is shown in the next theorem.

Theorem 4.1**.**

For sufficiently large $L$ there exists a constant $q<1$ such that

[TABLE]

for all $n$ and, hence, the error at the $n$-th iteration is

[TABLE]

Analogous statements are true for $\Psi_{+}^{(0)n}$.

Proof.

First consider (30) on the line $x+ai$. Let $\max_{x\in\mathbb{R}}|(K_{1})_{-}^{-1}(x+ai)|=c_{1}$, note that $K_{1}$ is non-zero on the strip so the constant is well defined. Then, using (32), we have

[TABLE]

This implies that

[TABLE]

Hence we obtained the desired result with $q=d_{1}d_{2}\epsilon_{2}\epsilon_{1}$ using (42). The contraction mapping theorem can be applied if $q<1$ to show that the iterations converge to the exact solution $\Phi_{-}^{(L)}$. Note that we can make $\epsilon_{1}$ and $\epsilon_{2}$ (and hence $q$) arbitrarily small by taking $L$ sufficiently large.

Similarly, we consider (31) on the line $x+bi$. In the same manner we derive

[TABLE]

∎

Note that the procedure converges for all possible initial functions. But the accuracy of the initial step determines how fast the small desired error will be achieved.

5. Examples

In this section the method proposed in this paper will be illustrated numerically. Two examples will be considered, in both cases the exact solutions are known and will be compared with the outcomes of the iterative procedure. The first example is of the type (1) and is chosen to be as simple as possible. The second example has been studied by other researchers in connection to integral equations. The resulting Wiener–Hopf system is more general than (1), but it can also be reduced to solving equations similar to (22) and (23) and hence the derivations in this paper apply. A more involved examples of application of the proposed method is found in [my+Lorna] and is used in the setting of acoustics, see Section 2.3 for more details. Note that in the latter case no exact solution is known.

Example 1**.**

The first numerical example will be the simplest possible in order to illustrate the theory. Consider (1) with

[TABLE]

where $\lambda$ is a complex parameter and $B_{+}(\alpha)$ is an arbitrary function satisfying the conditions stated in Section 2.1. For the forcing terms take

[TABLE]

In this example $a=b=\Re(\lambda)$ and so we require $\Re(\lambda)>0$. We have that

[TABLE]

Then equation (22) and (23) become

[TABLE]

and we are required to find $\Phi_{-}^{(L)}(\alpha)$ and $\Psi_{+}^{(0)}(\alpha)$. In this simple case we can exactly solve these coupled equations and find constants $k_{1}$ and $k_{2}$ explicitly. Define

[TABLE]

then the constants are given by

[TABLE]

Next, we will solve these coupled Wiener–Hopf equations using an iterative procedure described in Section 3. The first step is to neglect the the term with $e^{-i\alpha L}$ in (46) and hence uncouple the equations and obtain a value for $k_{2}^{(0)}$. This value is then substituted into (47) to obtain $k_{1}^{(0)}$. The iterative procedure leads to the following values

[TABLE]

This converges to the actual solution as long as $b<1$. In most cases the convergence is very fast and the line of the first iteration is indistinguishable to the approximate solution. In the cases when the convergence is slow (small $L$ and small $\Re(\lambda)$) the first iteration still retains some features of the solution, see Figure 1 and 2. The first guess has the correct overall shape and all the consecutive iterations have the maxima and the minima in the right places even when the magnitude is quite different. Also note that even in the cases of slow convergence, the second iteration is already close to the actual solution.

In this special case it is possible to say more about the convergence of the solution compared to the general case (Section 4). It is easy to see that

[TABLE]

Hence

[TABLE]

This means that

[TABLE]

This can be verified numerically and is illustrated in Figure 3.

Example 2**.**

The next example is more complicated and arises from an integral equation. The ability to solve integral equations is important in many applications [Vekua, Musk]. Consider the following (one-sided) integral equation

[TABLE]

with the matrix kernel given by

[TABLE]

with $f(x)$ a forcing function, $\lambda$, $L$ real parameters and $u(x)=(u^{(1)}(x),u^{(2)}(x))^{T}$ to be determined. This system was considered in [Abr_exp] and more recently in [Antipov]. This integral equation can be reduced to a Wiener–Hopf equation by extending the range of $x$ to the whole real line and applying the Fourier transform. The resulting Wiener–Hopf equation is of a more general type than (1). It has been shown [Antipov] that the solution could be reduced to finding two constants $C_{1}$ and $C_{2}$ by the Wiener–Hopf method. In this example it is possible to obtain the exact solution, which provides a good way to test out the ideas that are introduced in this paper. Once the solution to the Wiener–Hopf equation is obtained, the inverse Fourier transform will provide the solution to the integral equation.

First we will need to define some functions. Let us take two constants $\lambda_{0}=\sqrt{1-2\lambda}$ and $\lambda_{1}=\sqrt{1-4\lambda}$ for some $\lambda\in(-\infty,0.25]$, then we define

[TABLE]

The forcing is taken as $F_{1}^{+}(\alpha)=\frac{1}{\alpha-2i}$ and $F_{2}^{+}=0$ and we use the following additive splittings

[TABLE]

The Wiener–Hopf method reduces the solution to equations similar to (22) and (23). The unknown functions $U_{-}^{(2)}(\alpha)$ and $U_{+}^{(1)}(\alpha)$ are the half-range Fourier transforms (8) of $u^{(2)}(\alpha)$ and $u^{(1)}(\alpha)$. They are given by

[TABLE]

The constants can be found explicitly

[TABLE]

where

[TABLE]

By employing the iterative procedure we obtain

[TABLE]

In fact, each iteration adds two more terms in the Taylor expansion of the constants. For example, for $\lambda=0.1$ and $L=0.0001$ the maximum error for $U_{-}^{(2)0}$ is $10^{-4}$ and for $U_{-}^{(2)1}$ is $10^{-6}$. Note that the convergence speed of iterations depends on $b^{2}$, in the same way as in the previous example. It fact, it is small for all values of $\lambda$ and $L$ and $b^{2}\leq 0.18$. This shows that in this example the convergence is extremely fast for all values of the parameters. Another factor, which influences the error of the iterative solution, is suitability of the initial value. For example, with $\lambda=-15$ and $L=0.04$ the initial guess $U_{+}^{(1)0}$ is very bad, see Figure 4 the red dotted line. But since $b^{2}=0.0747$, even the first iteration $U_{+}^{(1)0}$ is already very accurate, see Figure 4 the blue dashed line.

6. Conclusion

The Wiener–Hopf method is a powerful tool for solving boundary value problems, which has been applied in an impressive array of situations. In the case of scalar Wiener–Hopf equations the solution is algorithmic. In the matrix case, the matrix Wiener–Hopf factorisation cannot be in general obtained and hinders the use of the method. The present paper presents an algorithmic way of bypassing this step, for a class of matrix functions (1). Only scalar Wiener–Hopf splittings are used in an iterative procedure. We provide the conditions for the convergence of iterations to the exact solution. This also enables us to estimate the error at each iteration. The numerical examples show that in most cases only a few iterations are required.

It is clear that this method could be applied to a wider class of Wiener–Hopf systems than (1). For instance, the matrix in Example 2 is not triangular, however it still can be solved using our methods. Further work could be done to extend this method to non-triangular matrix functions with exponential factors.

This method has been motivated by applications and was already used in [my+Lorna]. There is scope for using this method in a variety of boundary value problems with finite geometries or more than one transition point in the boundary conditions. Possible applications could be in different areas such as electromagnetism, fracture mechanics and economics. It has also been shown in numerical Example 2 that some integral equations can also be solved, this even opens a larger scope of applications.

7. Acknowledgements

I am grateful to Professor Peake for suggesting and discussing the acoustic problem which motivated this research. I also benefited from useful discussions with Professors Mishuris and Rogosin. I acknowledge support from Sultan Qaboos Research Fellowship at Corpus Christi College at University of Cambridge.