# Computing the oscillatory Cauchy operator¶

In [9]:
using ApproxFun, ApproxFunRational, ApproxFunOrthogonalPolynomials, ApproxFunFourier, Plots
# In the Julia REPL type "]" followed by "dev https://github.com/tomtrogdon/ApproxFunRational.jl"
# to install ApproxFunRational
# This will eventually be unified with RiemannHilbert.jl and SingularIntegralEquations.jl


In applications, many singular integral equations on $\mathbb R$ are of the form

$$$$u(s) - \begin{bmatrix} r_{11}(s)e^{i \alpha_{11} s }& r_{11}(s)e^{i \alpha_{12} s }\\ r_{21}(s)e^{i \alpha_{21} s } & r_{22}(s)e^{i \alpha_{22} s } \end{bmatrix} \mathcal C_{\mathbb R}^- u(s) = \begin{bmatrix} b_{1}(s)e^{i \alpha_{1} s }\\ b_{2}(s)e^{i \alpha_{2} s } \end{bmatrix}$$$$

where $u : \mathbb R \to \mathbb C^{1 \times 2}$.

• These usually arise as the singular integral reformulation of a Riemann-Hilbert or Weiner-Hopf problem.

• If all functions $r_{ij}$ and $b_j$ decay rapidly and the constants $\alpha_{ij}, \alpha_j$ are not too large, the methods we have developed so far can work.

• We will soon encounter situations where these restrictions do not hold and we then need to develop new methods to solve such an integral equation. One approach uses a so-called "oscillatory Cauchy transform".

The package ApproxFunFourier.jl includes functionality to perform an expansion of a function in the basis

$$$$f(x) \approx \sum_{j=-N}^N c_j \left( \frac{ i + x}{i-x} \right)^j.$$$$

In [10]:
j = -3; f = x -> ((1im + x)/(1im - x))^j
F = Fun(f,Laurent(PeriodicLine()))
F.coefficients

Out[10]:
7-element Array{Complex{Float64},1}:
-9.71445146547012e-17 - 1.2969353840194616e-16im
-3.3643862022641906e-16 + 1.2218701970925527e-18im
-2.3443003056032104e-17 + 3.60292215344155e-18im
-8.341116186757556e-17 + 1.542481527072819e-16im
-1.741099146248869e-16 + 5.849842599589495e-17im
0.0 + 4.981714270155072e-17im
1.0 - 4.4254798518585714e-16im
In [11]:
j = -1; f = x -> ((1im + x)/(1im - x))^j
F = Fun(f,Laurent(PeriodicLine()))
F.coefficients

Out[11]:
3-element Array{Complex{Float64},1}:
-5.551115123125783e-17 + 1.4592936398578187e-17im
0.0 - 7.15148590763743e-19im
1.0 - 1.4090699789001261e-16im

This expansion is computed by sampling $f$ at a set of points $x_0,x_1,\ldots,x_{2N}$ (evenly spaced on $[0,2\pi)$ and then mapped via an arctan transformation) and applying the fast Fourier transform to compute the coefficients.

Given $$$$f(x) \approx \sum_{j=-N}^N c_j \left( \frac{ i + x}{i-x} \right)^j.$$$$ we can define the function $$$$e^{i \alpha x} f(x) \approx \sum_{j=-N}^N c_j e^{i \alpha x} \left( \frac{ i + x}{i-x} \right)^j.$$$$ using ApproxFunRational.jl

In [12]:
f = x -> 1/(1 + x^2 + x^4); α = 2.
F = Fun(zai(f),OscLaurent(α)) # zai = zero at infinity
x = 1.; f(x)*exp(1im*α*x) - F(x)

Out[12]:
2.7755575615628914e-17 - 1.6653345369377348e-16im

If we want to, instead, use the basis $$$$f(x) \approx \sum_{j=-N}^N c_j \left( \frac{ Li + x}{Li-x} \right)^j.$$$$ for $L > 0$, we can:

In [13]:
j = -1; L = 2.; f = x -> ((1im*L+ x)/(1im*L - x))^j
F = Fun(f,OscLaurent(0.0,L)) # Or Laurent(PeriodicLine{false,Float64}(0.0.,L))
F.coefficients

Out[13]:
4-element Array{Complex{Float64},1}:
-2.7755575615628907e-17 + 3.061616997868383e-17im
0.0 + 5.265055686820291e-17im
1.0 - 1.1388289682557058e-16im
2.7755575615628907e-17 + 3.061616997868383e-17im 

Now, if $f(\infty) = 0$ then we enforce $$$$\sum_{j=-N}^N c_j ( -1)^j = 0$$$$ and we can write $$$$f(x) \approx \sum_{j=-N}^N c_j \left [\left( \frac{ Li + x}{Li-x} \right)^j - (-1)^j \right].$$$$ And to keep consistent with other notation we write $$$$f(x) \approx \sum_{j=-N}^N c_j (-1)^j \left [\left( \frac{ x + Li }{x - Li} \right)^j - 1 \right].$$$$

### The oscillatory Cauchy operator¶

The oscillatory Cauchy operator $\mathcal C_\alpha^{\pm}$ is defined to be $$$$\mathcal C_\alpha^{\pm} f(x) = \lim_{\epsilon \downarrow 0}\frac{1}{2 \pi i} \int_{\mathbb R} e^{i \alpha x'} f(x') \frac{dx'}{x' - (x \pm i \epsilon)}, \quad \mathcal C^\pm := \mathcal C_0^\pm.$$$$

This leads naturally to the question of how to compute $$$$\lim_{\epsilon \downarrow 0}\frac{1}{2 \pi i} \int_{\mathbb R} \underbrace{e^{i \alpha x'} \left[\left( \frac{ x' + Li }{x' - Li} \right)^j - 1 \right]}_{R_{-j,\alpha}(x'/L)} \frac{dx'}{x' - (x \pm i \epsilon)}.$$$$ Note that $R_{j,\alpha}(x'/L)$ has a pole at $x' = i L\mathrm{sign}(j)$.

This is interesting question and it will lead us right back to orthogonal polynomials.

From [TT 2013]: If $\alpha j \leq 0$ then \begin{align*} \mathcal C^+ R_{j,\alpha}(x) &= \begin{cases} R_{j,\alpha}(x) & j > 0, \\ 0 & j < 0, \end{cases}\\ \mathcal C^- R_{j,\alpha}(x) &= \begin{cases} 0 & j > 0 \\ -R_{j,\alpha}(x) & j < 0.\end{cases} \end{align*} There exists a function $\eta_{j,n}(\alpha)$, $\alpha j > 0$ such that \begin{align*} \mathcal C^+ R_{j,\alpha}(x) & = \begin{cases} - \displaystyle \sum_{n=1}^{j} \eta_{j,n}(\alpha) R_{n,0}(x) & j > 0,\\ R_{j,\alpha}(x) + \displaystyle\sum_{n=1}^{j} \eta_{j,n}(\alpha) R_{-n,0}(x) & j < 0, \end{cases}\\ \mathcal C^- R_{j,\alpha}(x) & = \begin{cases} -R_{j,\alpha}(x)- \displaystyle\sum_{n=1}^{j} \eta_{j,n}(\alpha) R_{n,0}(x) & j > 0,\\ \displaystyle\sum_{n=1}^{j} \eta_{j,n}(\alpha) R_{-n,0}(x) & j < 0. \end{cases} \end{align*}

As a first attempt, the coefficients $\eta_{j,n}(\alpha)$ were written as a sum of hypergeometric functions after a pure residue calculation. But after some thought, one realizes that, for example, $$\displaystyle \sum_{n=1}^{j} \eta_{j,n}(\alpha) R_{n,0}(x)$$ is a non-oscillatory function and the coefficients can be computed if the function itself can be computed on the appropriate grid $x_0,x_1,\ldots,x_{2N}$ (to then apply the fast Fourier transform).

\begin{align*} \sum_{n=1}^{j} \eta_{j,n}(\alpha) R_{\sigma n,0}(x) &= r_{j,\alpha}\left( \frac{-2 i \sigma L }{x + \sigma i L} \right) :=\mathrm{Res}\,\left\{ R_{j,\alpha}(x'/L) \frac{1}{x' - x} ; x' = -\sigma iL \right\} \\ &= \sum_{n = 1}^{|j|} \gamma_{j,n}(\alpha) \left( \frac{-2 i \sigma L }{x + \sigma i L} \right)^n, \sigma = \mathrm{sign}(j), \end{align*} where \begin{align*} \gamma_{j,n}(\alpha) = - e^{- |\alpha| L} L_{|j|-n}^{(n)}( 2 |\alpha| L) \end{align*} and $L_{k}^{(\alpha)}(x)$ is the $k$th-order generalized Laguerre polynomial. We then use the well-known relation for Laguerre polynomials \begin{align*} L_n^{(\alpha)}(x) = L_n^{(\alpha+1)}(x) - L_{n-1}^{(\alpha + 1)}(x) \end{align*} to obtain the recurrence relation for $j \geq 1$ \begin{align*} r_{j,\alpha}(z) = (1 + z) r_{j-1,\alpha}(z) + z e^{- |\alpha| L}( L_{j-1}^{(1)}(2 |\alpha| L) - L_{j-2}^{(1)}(2 |\alpha| L)), \quad z = \frac{-2 i \sigma L }{x + \sigma i L}. \end{align*}

So, then the seemingly complicated sum takes the form: \begin{align*} \sum_{j=1}^N c_j \sum_{n=1}^{j} \eta_{j,n}(\alpha) R_{\sigma n,0}(x) = \sum_{j=1}^N c_j r_{j,\alpha}(z) \end{align*} and if this is computed on the grid $x_0,x_1,\ldots,x_{2N}$ ( $O(N^2)$ flops because of the recurrence ) the fast Fourier transform will give the coefficients $d_j(\alpha)$ in the expansion $$\sum_{j=1}^N c_j r_{j,\alpha}(z) = \sum_{j = 1}^N d_j(\alpha) R_{\sigma j,\alpha}(x)$$

In [17]:
f = x -> 1/(x + 2im) + exp(-x^2)
F = Fun(zai(f),OscLaurent(2.0,1.0))
𝓒⁺ = Cauchy(+1)
𝓒⁻ = Cauchy(-1)
@time (𝓒⁺*F)(.1)
@time (𝓒⁺*F)(.1) - (𝓒⁻*F)(.1) - F(.1)

  0.006612 seconds (18.33 k allocations: 1.279 MiB)
0.008125 seconds (18.50 k allocations: 1.695 MiB)

Out[17]:
0.0 + 0.0im

With sufficient decay, the relation \begin{align*} \lim_{z \to \pm i \infty} z\int_{\mathbb R} f(x') \frac{d x'}{x' -z} = -\int_{\mathbb R} f(x')d x' \end{align*} gives us the Fourier transform of $R_{j,0}(x)$: \begin{align*} \int_{\mathbb R} e^{- i \alpha x'} R_{j,0}(x') {d x'} = \begin{cases} - 4 \pi L e^{-|\alpha| L} L_{|j|-1}^{(1)}( 2 |\alpha| \sigma) & j \alpha > 0, \\ -2 \pi |j| L & \alpha = 0,\\ 0 & j \alpha < 0, \end{cases} \end{align*} and a numerical method for computing Fourier transforms.

Note: From [Webb & Iserles 2019] this implies that the basis $\{R_{j,0}\}$ has a skew-symmetric differentation matrix!

In [18]:
f = x -> 1/(x - 2im) + exp(-x^2)
F = Fun(zai(f),OscLaurent(2.0,1.0))
𝓕 = FourierTransform(-1.0)
(𝓕*F)(.1)

Out[18]:
0.5885254797806003 + 0.09421998795110886im
In [19]:
F̂ = 𝓕*F  # This uses Laguerre functionality in ApproxFunOrthogonalPolynomials.jl
ω = -5:.01:5
y = map(F̂,ω)
plot(ω,real(y))
plot!(ω,imag(y))

Out[19]: