## Abstract

This paper presents a novel method of approximating the scalar Wiener–Hopf equation, and therefore constructing an approximate solution. The advantages of this method over the existing methods are reliability and explicit error bounds. Additionally, the degrees of the polynomials in the rational approximation are considerably smaller than in other approaches. The need for a numerical solution is motivated by difficulties in computation of the exact solution. The approximation developed in this paper is with a view of generalization to matrix Wiener–Hopf problems for which the exact solution, in general, is not known. The first part of the paper develops error bounds in *L*_{p} for . These indicate how accurately the solution is approximated in terms of how accurately the equation is approximated. The second part of the paper describes the approach of approximately solving the Wiener–Hopf equation that employs the rational Carathéodory–Fejér approximation. The method is adapted by constructing a mapping of the real line to the unit interval. Numerical examples to demonstrate the use of the proposed technique are included (performed on Chebfun), yielding errors as small as 10^{−12} on the whole real line.

## 1. Introduction

The Wiener–Hopf method is used for a broad collection of PDEs that arise in acoustic, finance, hydrodynamic, elasticity, potential and electromagnetic theories [1,2]. It is an elegant method based on the exploitation of the analyticity properties of the functions. For a scalar Wiener–Hopf equation, the solution can be expressed in terms of a Cauchy-type integral [3, ch. 1.3].

In more complicated scalar Wiener–Hopf equations, the exact solution is difficult or slow to compute [4–6]. Approximate solutions were considered early on, but they were mainly constructed using ad hoc observations [3, ch. 4.5]. In 2000, a systematic method of approximating the Wiener–Hopf equations was published by Abrahams [7]. Since then, it has proved popular and found applications in different branches of mathematics, including finance [8].

The method proposed in [7] is based on uniform approximations of the kernel on the whole strip by a two-point Padé approximation (the two points being 0 and on the real axis). However, there are two issues that make the application of this method difficult. The first is that it is unclear when the two-point Padé approximation has small errors on the whole strip [9]. Secondly, even if the maximum error of approximating the kernel on the whole strip is known, there are no error bounds presented to calculate the resulting error in the solution.

Another motivation for the development of approximate methods is the matrix Wiener–Hopf problem. The determination of a good numerical solution is important for the matrix Wiener–Hopf problem, as there is not yet a constructive method of solving it in general.

This paper aims to present a consistent method of approximately solving the scalar Wiener–Hopf equation. The following questions are addressed.

### Problem 1.1

Given a scalar Wiener–Hopf problem, find an approximate solution that can be demonstrated to be within a given accuracy from the exact solution.

### Problem 1.2

Perform problem 1.1 computationally in a way that is reliable, optimal and easy.

Problem 1.1 will be addressed in §3 with new estimations for the error in the Wiener–Hopf factors. This is done though expressing the factors in terms of the Hilbert transform, which is connected to the Cauchy singular integral. Note that although the method proposed in this paper involves the construction of an approximate solution to the given Wiener–Hopf equation, the solution is the exact solution of a perturbed Wiener–Hopf equation.

The main difficulty in problem 1.2 is that the Wiener–Hopf equation is set on an unbounded interval. We propose the novel method based on an appropriate mapping to the unit interval, see §4. This allows the Carathéodory–Fejér rational approximation to then be used in §5, and numerical examples are given in §6.

## 2. Preliminaries

The following conventions will be used throughout the paper. *A strip* around the real axis (given *τ*_{−}<0<*τ*_{+}) is . The subscript + (or −) indicates that the function is *analytic in the half-plane τ*>*τ*_{−} (or *τ*<*τ*_{+}). Functions in the Wiener–Hopf equation without the subscript are analytic in the strip. The Wiener–Hopf problems are recalled below.

The *multiplicative Wiener–Hopf* problem is: given a function *K* (analytic, zero-free and *K*(*α*)→1 as in the strip), to find functions *K*_{+} and *K*_{−} that satisfy the following equation in the strip:
2.1In addition, *K*_{−} and *K*_{+} are required to be:

— analytic and non-zero in the respective half-plane and

— of

*subexponential growth*of the respective half-planes,

The function *K*(*α*) in the multiplicative Wiener–Hopf factorization will be referred to as *the kernel*. The functions *K*_{+} and *K*_{−} will be called *factors* of *K*. These factors are unique up to a constant. In other words, if there are two such factorizations, *K*=*K*_{+}*K*_{−} and *K*=*P*_{+}*P*_{−}, then
where *c* is an arbitrary complex number [10]. In this paper, they will be normalized so that *K*_{±}(*α*)→1 as in the strip.

### Remark 2.1

The above normalization at infinity allows control of the regularity of the Fourier transform of an approximation. This is useful if a factorization of the Fourier transform is used, say, for a solution of differential equations.

The multiplicative Wiener–Hopf factorization can be reduced to the *additive Wiener–Hopf* problem via application of the logarithm. Under the conditions on the kernel *K*, let , then the additive Wiener–Hopf problem is: given a function *f*, find functions *f*_{+} and *f*_{−} that satisfy the following equation in the strip:
2.2Additive and multiplicative problems are equivalent only in the scalar case. In the matrix Wiener–Hopf problem, *K*_{+} and *K*_{−} will in general not be commutative, so no such equivalence is possible.

The existence of solutions to the additive and multiplicative scalar Wiener–Hopf problem is addressed in, for example, Noble [3, ch. 1.3].

These two Wiener–Hopf splittings are key to solving the *general Wiener–Hopf problem*. That is, given functions *A* and *C*, find functions *Φ*_{+} and *Ψ*_{−} that satisfy the following equation in the strip:
2.3For more details about the Wiener–Hopf problem, see [3].

### (a) The Riemann–Hilbert problem

The Wiener–Hopf problem is a special case of the Riemann–Hilbert problem. Roughly, the Riemann–Hilbert problem has a more general contour instead of the strip, and the conditions on the function *A*(*α*) are weakened. In particular, it is required that *A*(*α*) is only Hölder continuous on the contour.

Let *G*(*t*),*g*(*t*) be Hölder continuous functions on a simple contour *Σ*. The *Riemann–Hilbert problem* is to find an analytic function *E*(*α*) that has values *E*_{−}(*t*),*E*_{+}(*t*) on the contour as the limit is taken from different sides of the contour and that satisfies the equation
The Wiener–Hopf equation (2.3) can be considered as a special case of the Riemann–Hilbert problem so a good approximation of the values of the function on the real line is sufficient. This is simpler than trying to construct an approximation that agrees well on the whole strip of analyticity.

## 3. Estimates on the approximate Wiener–Hopf factorization

The approximate solution to (2.1) can be found by approximating the kernel in the Wiener–Hopf equation. This section quantifies the difference between the modified and the original Wiener–Hopf equation, to address problem 1.1 introduced earlier. More precisely, let , the aim is to bound in terms of *ϵ*_{p} and computable quantities of *K*(*α*).

The first step is to link the Wiener–Hopf factorization to the Hilbert transform *H*(*f*)(*y*) of *f* with . If
then [11, ch. 2]
and hence
The advantage of expressing *f*_{±} in terms of the Hilbert transform is the following theorem.

### Theorem 3.1 (Titchmarsh–Riesz [11, p. 92])

*Let* *for some* *then
**where the best constant is*
3.1

Based on these classical estimations, we obtain the following bounds on the additive Wiener–Hopf factorization.

### Lemma 3.2 (Additive bounds in *L*_{p} for )

*Let f*(*y*)=*f*_{+}(*y*)+*f*_{−}(*y*) *and* *with* , *then*
*where c*(*p*) *is defined as in* (*3.1*).

### Proof.

Expressing in terms of the Hilbert transform and using Minkowski’s inequality,
Here, the linearity property of the Hilbert transform, *H*(*f*)+*H*(*g*)=*H*(*f*+*g*) is used. □

The theorem gives bounds on the error on the real line, but, in fact, since the functions are analytic, the bounds will hold in the upper/lower half-planes by the maximum modulus principle.

It will be useful to express the Wiener–Hopf factors in terms of the Hilbert transform, The main new result of this section is the following.

### Theorem 3.3 (Multiplicative bounds in *L*_{p} for )

*Let K(y) and* *be two kernels and m<∥K∥*_{p}*<M. If* *, then*
*where c(p) is defined in (3.1).*

### Proof.

The first step is to express the product factorization *K*(*y*)=*K*_{+}(*y*)*K*_{−}(*y*) as an additive factorization. Taking logarithms,
From the assumptions on the kernel *K*(*y*), will be a continuous and single-valued function, and hence the additive decomposition is well defined. Given
and then using the mean value inequality (and ),
Applying lemma 3.2 gives
Finally, note that as
and so
from which the result follows. □

### Remark 3.4 (Real kernels)

In the case when kernels *K* and are real valued, the bounds could be simplified, since then and so

### Remark 3.5 ( norm)

The above theorem does not include the case of the norm, in fact, no such result is possible for the norm, as is demonstrated by the counter example below.

Consider *K*_{1}(*y*)=1 defined on the real line and
and
Then,
but it has been shown by making *n*_{2}−*n*_{1} small that the factors can differ by an arbitrary large amount in the norm [12].

Below, the proposed method of the solution of problem 1.1 is presented.

— Approximate, with arbitrary accuracy, the Wiener–Hopf kernel

*K*by rational functions , which can be done by Runge’s approximation theorem [13].— Perform the Wiener–Hopf factorization of the rational kernel by inspection. The will have the form where the

*z*_{i}(*p*_{j}) are all the zeros (poles) of that lie in the lower half-plane. The other factor will have the remainder poles and zeros in the upper half-plane.— Using theorem 3.3, the error can be calculated.

The rest of the paper will concentrate on problem 1.2, i.e. the practical aspect of problem 1.1.

## 4. Mapping of the real line

In §3, it was proved that the approximation of the kernel results in a computable error in the Wiener–Hopf factors. Thus, to simplify the problem, one may wish to approximate a kernel by another that is easier to factorize, e.g. by rational functions. To construct such an approximation, the first step is to employ a mapping of the real line to the unit circle or the interval [−1,1]. Those mappings transform the problem of approximation on the whole real line to a simpler problem on the unit circle or the interval. Methods for approximating on the unit circle and the interval will be discussed in §5.

To easily distinguish between functions defined on the real axis, the unit circle and the interval [−1,1], the following notation will be used:

— functions on the real line will be denoted by a capital letter and variable

*y*, for example,*F*(*y*);— functions on the unit circle will be denoted by a lower case letter and variable

*z*, for example,*f*(*z*);— functions on the interval [−1,1] will be denoted by bold capital letters and variable

*x*, for example,**F**(*x*); and— the rational approximation for functions will be denoted by letters with a tilde, for example,

It is typical in applications for the kernel to be an even function (i.e. *F*(*y*)=*F*(−*y*)), and the mapping is simplified in this case. Therefore, it is instructive to consider the mapping for even functions before considering the general case.

### (a) Even functions

Define a conformal Möbious map *M*(*y*) that maps the real line to the unit circle as follows [14, ch. 13]:
Explicitly, it is given by
The next map is a projection of the unit circle to the interval [−1,1] by the Joukowski map (a conformal map that, incidentally, comes from aerodynamics [15, ch. 5.2]). The *Joukowski map* is
4.1Restricted to the unit circle, the map returns the real part of *z*. This step requires the function to be even. Even functions will be mapped to a function on the unit circle with the property *f*(*z*)=*f*(*z*^{−1}). Composing these maps together gives
This enables the construction of a function **F**(*x*) on the interval from a given function *F*(*y*) on the real line as follows:
Then, the algorithm of §5 can be used to rationally approximate **F**(*x*) on [−1,1] by . Suppose that the error in this approximation in is *λ*.

Then, map the function back to the real line by
Furthermore, is a rational function on the real line that approximates *F*(*y*) on the real line with the error at most *λ*.

The maps described above are combined with the Carathéodory–Fejér algorithm in Chebfun to give the code used for the numerical examples in §5.

### Remark 4.1

If the kernel is odd, i.e. *F*(*y*)=−*F*(−*y*), the kernel can be squared to produce an even kernel. Then, the approximation can be produced as above and square rooted. Note that the explicit multiplicative Wiener–Hopf factorization of the square root of a rational function is carried out by inspection in the same way as for the rational function. Also note that any function can be written as a sum of even and odd functions, thus the above extends to a general function.

### (b) General functions

For a general kernel, the above map can be modified. The above method contains the following ideas. The real line is conformally mapped to the circle. Then, only half of the circle is mapped to the interval [−1,1]; however, this poses no problem since the other half is the same owing to the function being even. The function is approximated and mapped back to the unit circle. The modification is needed to ensure that half of the function is not discarded. Thus, after mapping the real line to the circle, apply the map *S*: *z*→*z*^{2}; this makes the values run twice as fast and the whole function now fits on the half circle. Next, as before, map this half circle to the interval, approximate and map back. The values that are taken on the half circle are spread out to the circle by the map . Calculating the resulting map gives (figure 1)
4.2

### Remark 4.2

Before using a map, it is beneficial to apply a rescaling and a shift (i.e. a Möbious map) to the real line so that the part of the function that has a fast changing gradient fits in the interval [−5,5]. This part of the real line is well resolved (figure 1).

### (c) Associated orthogonal rational functions

In this section, the properties of the map (4.2) will be studied. This will be carried out by considering the new basis functions that are created by this map. A starting point is to choose a basis function for the circle; a good choice for this is the Fourier basis. Then, the Joukowski map will create the new basis for the interval [−1,1] as the image of the Fourier basis, under the change of coordinates [16], §1.6. This new basis consists of Chebyshev polynomials of the first kind (they are also the Fabe polynomials on [−1,1] up to a multiplicative factor [17]). Furthermore, the mapping *M*^{−1}*S*^{−1}*J*^{−1}(*x*) from [−1,1] to will create a new basis (from Chebyshev polynomials) called *rational Chebyshev functions*, *TB*_{n} [18]. They are defined as
where *T*_{n} are the Chebyshev polynomials of the first kind and is the map *JSM*(*y*) (see (4.2)). Note that, despite the name, *TB*_{n} are not all rational: *TB*_{2n+1} are rational functions divided by (*TB*_{2n} are all rational functions).

A useful property of *TB*_{n} is the shape of the domain of convergence. This is especially informative when the kernel has no singularities on the whole real line, including infinity. For standard basis functions, the domain of convergence has a specific shape, and the size is determined by the position of singularities [18]. For example, the Taylor series converges on a circle and the Chebyshev series converges on an ellipse, the sizes of which are dictated by the position of the singularities. The shape of the domain for *TB*_{n} is the exterior of two circles that are the image of lines parallel to the real axis under the Möbious map *i*(1+*z*)/(−1+*z*). This is derived from the corresponding domain of convergence for Chebyshev polynomials that, in turn, comes from the Fourier series [18].

Therefore, the base functions are well suited for functions that have singularities in the bounded part of the complex plane (as in example 6.1). This compares favourably to the domain of convergence of Taylor expansion for a function similar to the one in example 6.1. In particular, this explains why all methods of approximation, which are based on the Taylor expansion, do not produce good results on the whole real line (e.g. the one-point Padé approximation).

For functions with a singularity at infinity, the following theorem by Boyd is applicable.

### Theorem 4.3 (Convergence of the Fourier series for an algebraically decaying or asymptotically constant function [18, §5])

*If a function u(y) is free from singularities on the real axis and has the inverse power expansion*
4.3*as* *and a similar series as* *then the coefficients of its representation as a Fourier series in the new coordinate t (where z=e*^{2πit} *and z=SM(y), see (4.2)),*
4.4
*will have exponential convergence in the sense that |a*_{n}*| and |b*_{n}*| decrease with n faster than any finite inverse power of n.*

*In the case of (4.3) converging, then (4.4) converges geometrically.*

*In the case of (4.3) being asymptotic but divergent (then u(y) will be singular at infinity but all the derivatives will be bounded there), then (4.4) converges subgeometrically.*

### Remark 4.4

The theorem covers precisely the class of functions that arises as kernels in the additive and multiplicative Wiener–Hopf factorization. In the additive case, the kernel must decrease faster than some power and in the multiplicative case, they are asymptotically 1.

The above theorem shows one of the strengths of the *TB*_{n} compared with other base functions on an unbounded interval, e.g. the Hermite and sinc functions. The latter only converge algebraically for algebraically decaying functions.

## 5. Rational approximation

The best rational approximation in is not practical to use [19]. This section describes a near best method of rational approximation named AAK (Adamjan–Arov–Kreĭn) [20] in the general case and rational Carathéodory–Fejér [21,22] in the case of real-valued functions. The AAK theory is briefly reviewed at the end. Unfortunately, it seems that there is no software that performs AAK approximation; only software for rational Carathéodory–Fejér has been developed, Chebfun (Matlab) [19].

### (a) Real rational approximation

A brief summary of the method from [23] is presented below. Let *F*(*x*) be a real continuous function on [−1,1]. For any natural number *M*, *F*(*x*) has a partial Chebyshev expansion,
where
and *T*_{k}(*x*) are Chebyshev polynomials of the first kind.

An application of the Joukowski map (see (4.1)) to *F*_{M}(*x*) produces an associated function on the unit circle
The above is true because of the intimate connection of the Chebyshev polynomials and Joukowski transform given by
Define
where we are seeking a [*m*,*n*] rational approximation. Next, approximate *f*^{+}(*z*) on the unit circle by an extended rational function of the form
5.1where the numerator is a bounded analytic function in |*z*|>1 and the denominator has no zeros in |*z*|<1. Call a class of functions of the above form . Note that this is not a typical class to approximate by, but it is the one for which a neat solution exists. After the approximation (5.1) is obtained, it can be truncated to get a rational approximation. To obtain an approximation in , consider a real symmetric Hankel matrix
Let *λ*_{i} be eigenvalues of *X* arranged in decreasing order by magnitude of absolute value. And let (*u*_{1},…,*u*_{M+n−m}) be an eigenvector for the *λ*_{n+1} eigenvalue. Then, the following theorem is true.

### Theorem 5.1 (Takagi)

*The analytic function f*^{+} *has a unique best approximation* *on the unit circle |z|=1 in* *given by*
*where b is the finite Blaschke product*
*The approximation error in the* *norm on the unit circle is |λ*_{n+1}*|.*

Once the approximation is obtained, it is mapped back from the unit circle to the unit interval by the inverse of the Joukowski transform. A subsequent truncation gives a near best rational approximation on the interval.

This is a fast and efficient way of computing rational approximations. Furthermore, it enables prediction of the error on the approximation, even before the approximation is computed by considering eigenvalues.

### (b) Adamjan–Arov–Kreĭn approximations

This section presents more general results to the ones covered in §5*a*. These apply to complex-valued functions on the unit circle and are taken from [20].

### Definition 5.2

Given a natural number *n*, define as a class of bounded functions on the unit circle that can be expressed as
where *r*(*z*) is a rational function that has no more than *n* poles all inside the unit circle and .

The AAK approximation solves the following problem.

Given a function and *n* a natural number, find the best approximation in the norm from functions in .

In other words, it is necessary to find and a function (if it exists) such that The solution to the above is presented in the next theorem.

### Theorem 5.3 ([20])

*Let* *then D*_{n}*(f)=s*_{n}*(M) where M is the Hankel matrix built out of Fourier coefficients of f(z) and s*_{n} *is the nth singular value of it. Moreover,*
*where {ξ,η} is a Schmidt pair for a Hankel matrix M and where*

The next section illustrates theory with examples.

## 6. Numerical examples

In this section, numerical examples are provided to illustrate the theory (performed in Chebfun, Matlab).

### Example 6.1

The first example (figure 2) is
6.1Take the finite branch cut from i to *ki* and from −*i* to −*ki*. This kernel is closely associated with the matrix kernel factorization from problems in acoustics and elasticity, and was studied by Abrahams [7]. The case *k*=2 will be considered in the numerical examples. The first approach is to approximate this given function on an interval, a process called domain truncation. On the surface, the results look promising (figure 3*a*); the error decreases as the degree of the numerator polynomial is increased. Nevertheless, there are problems, the most obvious being that [20,4] or any other approximations plotted in figure 3*a* have very different behaviours at infinity than *F*(*y*). The error will be small, but only on the given interval, and will not be controlled outside it. It might be tempting to try to rectify this, since *F*(*y*)−1 will be zero outside a sufficiently large interval. One might suggest to look at [*n*,*m*] approximations where *n* is smaller than *m*: this type of function will go to zero at infinity. But, although the error on most of the real axis is very small, just outside the interval, it tends to be greater (and much larger than it is inside the interval). Another feature of this type of approximation is how the poles and zeros are positioned. Although some of the poles are positioned on the branch cuts [*i*,2*i*] and [−*i*,−2*i*], there are many singularities introduced elsewhere (figure 3*b*).

In example 6.2, the mapping constructed in §4 will be used to overcome the problems described. This will also result in smaller degrees of the polynomials.

### Example 6.2

The previous example 6.1 is treated with the mapping proposed in §4 for even functions. The function mapped to the interval [−1,1] becomes
In figure 4, the new function is plotted; note that in this form, it is much easier to approximate, since it has a slow changing gradient. This is confirmed by the error curve for the [4,4] approximation; the error is of order 10^{−9}. Furthermore, a better result is obtained with the [5,5] approximation with an error of order 10^{−12} (figure 5). At this point, the approximated function can be mapped back to the real line and the norm of the error will stay the same (see solid curves of figure 6). Note that it now becomes a [10,10] rational approximation. This example was chosen because the Wiener–Hopf factors can be easily seen by inspection. They are
The factors of the rational approximation can be easily calculated and compared with the exact solution above; this is plotted in figure 6 (dotted curve). This error is smaller than the error obtained in [7], even with [30,30] for the same function.

What is more, the poles and zeros also lie almost exactly on the branch cuts [*i*,2*i*] and [−*i*,−2*i*]. The close proximity of zeros and poles mimic the behaviour of a branch cut. This indicates that the behaviour in the whole complex plane is correct (figure 7). Note that the scale of the *x*-axis is 10^{−7}.

Bounds developed for real-valued kernels in §3 could be applied to this example. The *L*_{2} bounds will be used. Note that *c*(2)=1, *M*=1 and ; this gives
The *L*_{2} norm is easy to estimate numerically in Chebfun by using the *norm* command. For the [10,10] rational approximation, these calculations give
which agree with the bounds above.

### Example 6.3

The function is (figure 8)

This kernel is related to the one considered in Koiter’s paper of 1954 [24], and is typical for electrostatic and slow flow problems. In the original article, rational approximations were considered by choosing the coefficients by hand, and an accuracy of 10^{−2} was achieved. Once again, the function becomes easier to approximate (on the whole real line) once it is mapped to the interval [−1,1] (figure 8). The mapped function is (general mapping is used)
This can then be approximated with errors of 10^{−8} (giving a graph similar to figure 5) with [12,12], even though the point at infinity is singular (it is an accumulation point of a sequence of poles). Interestingly, the position of the zeros and poles (figure 9) is very similar to that found in [7], fig. 8 with a different method of approximating. This kernel has exact factorization that is given by
But, Matlab does not have an inbuilt complex *Γ* function, and hence the results of this approximation cannot be compared with the exact factorization. Even if there was an inbuilt complex *Γ* function to compare it with the approximate solution, the inbuilt function would need to have very high accuracy. This demonstrates that if the numerical values of the solution are needed, the approximate solution is as good as the exact solution.

## 7. Conclusion

This paper demonstrates, with the use of examples, the method of approximate solution of the Wiener–Hopf equation. The bounds presented allow prediction of the error in the approximation. The degrees of polynomials in the rational approximation are small, which allows simplification of the initial problem significantly. The next project would be to generalize the approach to the matrix Wiener–Hopf problem.

## Acknowledgements

I am grateful to Prof. Nigel Peake for suggesting this project and for support along the way. Anonymous referees made several useful comments that helped to improve this paper.

- Received December 10, 2012.
- Accepted February 22, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.