{'fraction_non_alphanumeric': 0.04438451341591626, 'fraction_numerical': 0.011244686200118847, 'mean_word_length': 4.245264924478542, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Constructive Method for Approximate Solution to Scalar Wiener-Hopf Equations

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 generalisation to matrix Wiener-Hopf for which the exact solution, in general, is not known. The first part of the paper develops error bounds in Lp for 1<p<\infty. These indicate how accurately the solution is approximated in terms of how accurate the equation is approximated. The second part of the paper describes the approach of approximately solving the Wiener-Hopf equation that employs the Rational Caratheodory-Fejer 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 error as small as 10^{-12} on the whole real line.


Introduction
The Wiener-Hopf method is used for a broad collection of PDEs which arise in acoustic, finance, hydrodynamic, elasticity, potential and electromagnetic theories [1,15]. It is an elegant method based on the exploitation of the analyticity properties of the functions. For the scalar Wiener-Hopf the solution can be expressed in terms of a Cauchy type integral [16,Ch. 1.3].
In more complicated scalar Wiener-Hopf equations the exact solution is difficult or slow to compute, see e.g. [2,14,24]. Approximate solutions were considered early on but they were mainly constructed using ad hoc observations [16,Ch. 4.5]. In 2000 a systematic way of approximating the Wiener-Hopf equations was published by I. D. Abrahams [3]. Since then it proved popular and found applications in different branches of mathematics including finance [10].
The method proposed in [3] 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 error on the whole strip [6]. 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.
Remark 2.1. The above normalisation at infinity allows to control the regularity of the Fourier transform of an approximation. This is useful if a factorisation of the Fourier transform is used, say, for a solution of differential equations.
The multiplicative Wiener-Hopf factorisation can be reduced to the additive Wiener-Hopf problem via application of logarithm. Under the conditions on the kernel K, let f = log K then additive Wiener-Hopf problem is: given a function f to find functions f + and f − which satisfy the following equation in the strip: (2.2) Additive and multiplicative problems are equivalent only in the scalar case. In the matrix Wiener-Hopf 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 [16,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 Ψ − , which satisfy the following equation in the strip: For more details about the Wiener-Hopf problem see [16].
2.1. 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(α) which has values E − (t), E + (t) on the contour as the limit is taken from different sides of the contour and which satisfies the equation: The Wiener-Hopf equation (2.3) can be considered as the 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.

Estimates on the Approximate Wiener-Hopf Factorisation
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 introduced earlier.
The first step is to link the Wiener-Hopf factorisation to the Hilbert Transform H(f )(y) of f with y ∈ R. If: then [17,Ch. 2]: and hence: The advantage of expressing f ± in terms of the Hilbert transform is the following theorem.
Theorem 3.1 (Titchmarsh-Riesz). [17, p. 92] Let f ∈ L p (R) for some 1 < p < ∞, then: where the best constant is: Based on these classical estimations we obtain the following bounds on the additive Wiener-Hopf factorisation.
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: Theorem 3.3 (Multiplicative Bounds in L p for 1 < p < ∞). Let K(y) andK(y) be two kernels and m < ||K|| p < M . If ||K(y) −K(y)|| p < ǫ p then: where c(p) is defined in (3.1).
Proof. The first step is to express the product factorisation K(y) = K + (y)K − (y) as an additive factorisation. Taking logarithms: log K(y) = log K + (y) + log K − (y).
Remark 3.4 (Real Kernels). In the case when kernels K andK are real valued, the bounds could be simplified, since then ||K ± (y)|| p < (M + ǫ p ) 1/2 and so: Remark 3.5 (L ∞ norm). The above theorem does not include the case of L ∞ norm, in fact no such result is possible for the L ∞ norm as is demonstrated by the counter example below.
Then: ||K 2 (y) − K 1 (y)|| ∞ = ǫ, but it has been shown by making n 2 −n 1 small the factors can differ by an arbitrary large amount in L ∞ norm, see [5].
Below, the proposed method of the solution of Problem 1 is presented: • Approximate, with arbitrary accuracy, the Wiener-Hopf kernel K by rational functionsK, that can be done by Runge's approximation theorem [18].
• Perform the Wiener-Hopf factorisation of the rational kernelK by inspection. TheK + will have the form: where the z i (p j ) are all the zeroes (poles) ofK which lie in the lower halfplane. The other factorK − will have the remainder poles and zeroes in the upper half-plane.
The rest of the paper will concentrate on Problem 2 i.e. the practical aspect of Problem 1.

Mapping of the Real Line
In the previous section 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 which is easier to factorise 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 the next section.
At this point it is important to fix notation. There will be functions defined on the real axis, the unit circle, and the interval [−1, 1]. To easily distinguish between them the following conventions will be used: • Functions on the real line will be denoted by capital letter and variable y, for example F (y); • Functions on the unit circle will be denoted by small 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); • The rational approximation for functions will be denoted by letters with a tilde, for exampleF (y).
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.
The easiest way to picture the map is to think of the plane plus infinity as a Riemann Sphere and the map as a composition of two rotations of the Riemann Sphere. 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 which, incidentally, comes from aerodynamics [20]). The Joukowski map is: Restricted 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: 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 the next section can be used to rationally approximate F(x) on [−1, 1] byF(x). Suppose that the error in this approximation in the L ∞ is λ.
Then map the function back to the real line by: FurthermoreF (y) is a rational function on the real line that approximates F (y) on the real line with the L ∞ 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 the next section.
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 factorisation of the square root of a rational function is done 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 function, thus the above extends to a general function.

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 due 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 S −1 : z → √ z. Calculating the resulting map gives: JSM (y) = y see Figure 1. This is derived from the corresponding orthogonality property of the Chebyshev polynomials. This property can be used to calculate the coefficients in the expansion of a function in terms of the Rational Chebyshev Functions.
Another useful property of T B 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 [8]. For example, Taylor series converge on a circle and Chebyshev series converge on an ellipse, the sizes of which are dictated by the position of the singularities. The shape of the domain for T B n is the exterior of two circles which 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 which in turn comes from the Fourier series [8].
So in particular those base functions are well suited for functions which have singularities in the bounded part of the complex plane, like Example 6.1. This compares favourably to the domain of convergence of Taylor expansion for a function like 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. one point Padé approximation).
For functions with a singularity at infinity the following theorem by Boyd is applicable.  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) 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 factorisation. 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 T B n compared to other base functions on an unbounded interval e.g. the Hermite and sinc functions. The latter only converge algebraically for algebraically decaying functions.

Rational Approximation
The best rational approximation in L ∞ is not practical to use, see [23]. This section describes a near best method of rational approximation named AAK (Adamjan-Arov-Kreȋn) [4] in the general case and Rational Carathéodory-Fejér [11,21] in the case of real valued functions. The AAK theory is briefly reviewed at the end.
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: where the numerator is a bounded analytic function in |z| > 1 and the denominator has no zeroes in |z| < 1. Call a class of functions of the above formR m,n . 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 inR m,n 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: The analytic function f + has a unique best approximatioñ r on the unit circle |z| = 1 inR m,n given by: where b is the finite Blaschke product b(z) = λ n+1 z M u i + · · · + u M+n−m z M+n−m u M+n−m + · · · + u 1 z M+n−m−1 .
And the approximation error in L ∞ 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 to predict the error on the approximation even before the approximation is computed by considering eigenvalues.

AAK Approximation.
This section presents more general results to the ones covered in the previous section. These apply to complex valued functions on the unit circle and are taken from [4].
Definition 1. Given a natural number n, define H n ∞ a class of bounded functions on the unit circle which can be expressed as: where r(z) is a rational function that has no more than n poles all inside the unit circle and h(z) ∈ H ∞ .
The AAK approximation solves the following problem: Given a function f (z) ∈ L ∞ and n natural number, find the best approximation in the L ∞ norm from functions in the H n ∞ . In other words it is required to find: and a function h(z) ∈ H n ∞ (if it exists) such that: The solution to the above is presented in the next theorem.
where {ξ, η} is a Schmidt pair for a Hankel matrix M and where: The next section illustrates theory with examples.

Numerical Examples
In this section the numerical examples are provided to illustrate the theory. All the numerical experiments are performed in Chebfun (MATLAB).  Figure 2. The function to be approximated F (y) = (y 2 +1) (y 2 +4) .
6.1. Example 1a. The first example (see Figure 2) is Take the finite branch cut from i to ki and from −i to −ki. This kernel is closely associated with the matrix kernel factorisation from problems in acoustics and elasticity and was studied by I. D Abrahams [3]. 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 (see Figure 3); 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 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 zeroes are positioned. Although some of the poles are positioned on the branch cuts [i, 2i] and [−i, −2i], there are a lot of singularities introduced elsewhere (see Figure 4).
In the next section, the mapping constructed in Section 4 will be used to overcome the described problems. This will also result in smaller degree of the polynomials.
6.2. Example 1b. The previous example (6.1) is treated with the mapping proposed in Section 4 for even functions. The function mapped to the interval [−1, 1] becomes: In Figure 5 the new function is plotted; notice 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 [5,5] approximation with error of order 10 −12 (see Figure 6). At this point the approximated function can be mapped back to the real line and the L ∞ norm of the error will stay the same, see Figure 7 (solid green curves). 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  x 10 −12 Figure 6. The error in approximating F(x) by [5,5]. are: The factors of the rational approximation can be easily calculated and compared to the exact solution above; this is plotted in Figure 7 (dotted red curve). This error is smaller than the error obtained in [3]    What is more, the poles and zeroes also lie almost exactly on the branch cuts [i, 2i] and [−i, −2i]. The close proximity of zeroes and poles mimic the behaviour of a branch cut. This indicates that the behaviour in the whole complex plane is correct (see Figure 8). Note that the scale of the x-axis is 10 −7 .
6.3. Example 2. The function is ( Figure 9): This kernel is related to the one considered in Koiter's paper of 1954 [12] and is typical for electrostatic and slow flow problems. In the original article, rational approximations were considered by choosing the coefficients by hand, 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 9. The mapped function is (the general mapping is used): This can then be approximated with errors of 10 −8 (giving the graph similar to Figure 6) 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 zeroes and poles (see Figure 10) is very similar to that found in [3, Fig 8] with a different method of approximating. This kernel has exact factorisation which is given by: K ± (y) = e −iπ/4 Γ(1/2 − iy/π) √ π Γ(1 − iy/π) (i ± α) 1/2 .  Figure 9. The initial function to be approximated F (y) = √ (y 2 +1) y tanh y (top) and the mapped function to be approximated F(x) = 1 x tanh x √ 1−x 2 (bottom).
small which allow to simplify the initial problem significantly. The next project would be to generalise 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 which helped to improve this paper.  Figure 10. The poles and zeroes in approximating F (y) by [12,12] the poles are filled in red circles whereas the zeroes are green circles.