**NB** - This notebook was originally made, some years ago, by Steven G. Johnson. I took the pre 1.0 Julia code and updated it so as to teach myself and others some useful concepts about generated functions in Julia. It is based on the idea of generated functions and staged programming of chapter 5.7 of Jeff Bezanson's PhD thesis on the Julia Language.This transcription is shared with the author's permission. The original code is in Appendix B of the same thesis.

β Miguel Raz, 2020

Input: kernel function $K(r)$. Output "first integral" of $K$ :

$\mathcal{K}_n(X) = \int_0^1 w^n K(wX) dw$for $X \in (0,\infty)$. Some small set of $n$ values, to be determined, will be supplied at compile time.

$K$ will be specified as a subtype of `AbstractKernel`

.

For some kernel types, $\mathcal{K}_n$ is known analytically. (See Table I in Homer's paper.)

$K$ may depend on some compile-time numeric parameters, specified as type parameters.

If the $\mathcal{K}_n$ integral is

*not*known analytically, we may have to compute it numerically. This will be done at*compile time*by staged functions, with the numeric integration results used to compute the coefficients of a Chebyshev polynomial fit, which can then be compiled into an efficient polynomial approximation of $\mathcal{K}_n$. However, to compute a Chebyshev approximation for a function defined on $(0,\infty)$, we will have to perform a coordinate transform from $X \to (0,1)$, and the type of coordinate transformation will depend on how fast $\mathcal{K}_n$ decays asymptotically as $X\to\infty$. This decay rate can be*specified via the type*of`K`

.

The $\mathcal{K}_n$ function will be parameterized by a `FirstIntegral{K,n}`

type parameterized by an `AbstractKernel`

type `K`

and an integer `n`

.

`julia> abstract type AbstractKernel end`

Any kernel ~ X^P for X βͺ S and ~ X^Q for X β« S

`julia> abstract type PowerLawScaling{P,Q,S} <: AbstractKernel end`

`julia> immutable FirstIntegral{K<:AbstractKernel,N} end`

```
immutable PowerLaw{p} <: PowerLawScaling{p,p,p} end # rα΅ power law
function (::FirstIntegral{PowerLaw{p},n})(X::Number) where {p,n}
return p >= 0 ? X^p / (1 + n + p) : inv(X^(-p) * (1 + n + p))
end
F = FirstIntegral{PowerLaw{-1}, 3}()
F(3.7)
```

`julia> @code_llvm F(3.7) # best to run this one yourself!`

In a general multi-physics BEM package, one might conceivably have a user-specified kernel $K$ for which the first integral $\mathcal{K}_n$ is *not* known analytically. Performing the integral numerically at runtime would be too expensive, however, especially since the integrand may have an integrable singularity.

Instead, we will perform the integral $\mathcal{K}_n$ at *compile-time* (or, at least, outside the innermost BEM loops) for various $X$ and use these values to compute a *Chebyshev interpolating polynomial* $C(\xi)$. This polynomial will then be used to generate an efficient $\mathcal{K}_n(X)$.

There are three tricky points:

$\mathcal{K}_n(X)$ will probably be singular as $X\to 0$, which means we can't fit it directly to a polynomial. (This is not a problem for how $\mathcal{K}_n$ is

*used*, because $\mathcal{K}_n$ is always used for integration over domains that do not include $X=0$.) This will be dealt with by requiring the user to specify the degree $p$ of the singularity as $X \to 0$: i.e. $\mathcal{K}_n(X) = O(X^p)$ for $X\to 0$. We will factor out this singularity from $\mathcal{K}_n$ and fit the remaining (non-singular) function to a polynomial.We need $\mathcal{K}_n(X)$ for $X \in (0,\infty)$, whereas polynomial interpolation requires a finite interval, typically $(-1,1)$. We will handle this by choosing a coordinate mapping $\xi(X) \in (-1,1)$. Because such a coordinate mapping is necessarily singular, however, it will screw up the convergence of polynomial interpolation if we choose the wrong degree of singularity β we want a mapping such that $\mathcal{K}_n(\xi(X))$ is nonsingular (e.g. a low-degree polynomial in $\xi$) as $X\to\infty$. To choose this, the user will specify a degree $q$ of the decay rate as $X\to\infty$, i.e. $\mathcal{K}_n(X) = O(X^q)$ for $X\to \infty$.

$X$ is dimensionful (it is a physical distance within one of the triangles or other geometric elements of the BEM basis). Mapping it to a dimensionless $\xi$ inevitably involves choosing a scale $s$ of $X$. This $s$ should be user-specified (e.g. it can be the median diameter of the BEM elements). That is, $\mathcal{K}_n(X) \sim X^p$ for $X \ll s$ and $\mathcal{K}_n(X) \sim X^q$ for $X \gg s$.

In summary, the user will specify a `PowerLawScaling{p,q,s}`

kernel type parameterized by `p`

, `q`

, and `s`

. They will also define a `(::PowerLawScaling{p,q,s})(r::Number)`

method that computes $K(r)$.

The polynomial fit will be performed as follows, assuming $p \le 0$ and $q \le 0$. First, we let $L_n(X) = \mathcal{K}_n(X) / (s^p + X^p)$, which eliminates the $x\to 0$ singularity while still having $L_n \sim X^q$ for $X \gg s$. Second, we let $X = (1-\xi)^{1/q} - 2^{1/q}$ [or equivalently $\xi = 1 - (x+2^{1/q})^q$], which maps $\xi \in (-1,1)$ to $X \in (0,\infty)$, and has the property that $X^q \approx 1-\xi$ as $\xi\to 1$, so the coordinate remapping produces a nice degree-1 polynomial. Finally, we fit $L_n(\xi(X)) = C(\xi)$ to a Chebyshev polynomial $C$, and compute $\mathcal{K}_n(X)$ via $\mathcal{K}_n(X) = (s^p + X^p) \times L_n(\xi(X))$.

The following routines compute the coefficients $c_n$ of a Chebyshev interpolating polynomial $C(x) = \sum_{n=0}^{N-1} c_n T_n(x)$ for a function $f(x)$ on $(-1,1)$, where $T_n(x) = \cos(n \cos^{-1}x)$ are the first-kind Chebyshev polynomials.

We compute these coefficients $c_n$ by first evaluating $f(x)$ at the Chebyshev points $\cos\left(\pi\frac{n+1/2}{N}\right)$ for $n=0,\ldots,N-1$, for which the Chebyshev sum is equivalent to a type-III discrete cosine transform (DCT-III), so that the coefficients $c_n$ are computed by a DCT-II. These are *not* quite the typical Chebyshev points, which correspond to a DCT-I: the difference is that the DCT-I corresponds to the closed interval $[-1,1]$, i.e. it includes the endpoints, whereas our function may involve terms that blow up at the endpoints (although the overall function should be okay) so we don't want to evaluate it there.

We also provide a function `evalcheb`

to evaluate $C(x)$ for any $x\in(-1,1)$ by a Clenshaw recurrence, and a macro version `@evalcheb`

(analogous to `Base.@horner`

) that generates a completely inlined version of this recurrence for the case where $c$ is fixed.

```
julia> #Pkg.add("FFTW");
julia> using FFTW
```

$N$ chebyshev points (order N) on the interval $(-1,1)$

`julia> chebx(N) = [cos(Ο*(n+0.5)/N) for n in 0:N-1]`

$N$ chebyshev coefficients for vector of $f(x)$ values on $chebx$ points $x$

```
julia> function chebcoef(f::AbstractVector)
a = FFTW.r2r(f, FFTW.REDFT10) / length(f)
a[1] /= 2
return a
end
```

Given a function $f$ and a tolerance, return enough Chebyshev coefficients to reconstruct $f$ to that tolerance on $(-1,1)$

```
julia> function chebcoef(f, tol=1e-13)
N = 10
local c
while true
x = chebx(N)
c = chebcoef(float[f(y) for y in x])
# look at last 3 coefs, since individual c's might be zero
if maximum(abs,c[end:end-2]) < tol * maximum(abs,c)
break
end
N *= 2
end
vβ = maximum(abs,c) * tol
return c[1:findlast(v -> abs(v) > tol, c)] # shrink to minimum length
end
julia> function chebcoef(f, tol=1e-13)
N = 10
local c
while true
x = chebx(N)
c = chebcoef(Float64[f(y) for y in x])
# look at last 3 coefs, since individual c's might be zero
if maximum(abs,c[end:end-2]) < tol * maximum(abs,c)
break
end
N *= 2
end
vβ = maximum(abs,c) * tol
return c[1:findlast(v -> abs(v) > tol, c)] # shrink to minimum length
end
```

Given $cheb$ coefficients $a$, evaluate them for $x$ in $(-1,1)$ by Clenshaw recurrence

```
julia> function evalcheb(x, a)
isempty(a) && throw(BoundsError())
-1 β€ x β€ 1 || throw(DomainError())
bβββ = bβββ = zero(x)
for k = length(a):-1:2
bβ = a[k] + 2x*bβββ - bβββ
bβββ = bβββ
bβββ = bβ
end
return a[1] + x*bβββ - bβββ
end
# inlined version of evalcheb given coefficents a, and x in (-1,1)
julia> macro evalcheb(x, a...)
isempty(a) && throw(BoundsError())
# Clenshaw recurrence, evaluated symbolically:
bβββ = bβββ = 0
for k = length(a):-1:2
bβ = esc(a[k])
if bβββ != 0
bβ = :(muladd(t2, $bβββ, $bβ))
end
if bβββ != 0
bβ = :($bβ - $bβββ)
end
bβββ = bβββ
bβββ = bβ
end
ex = esc(a[1])
if bβββ != 0
ex = :(muladd(t, $bβββ, $ex))
end
if bβββ != 0
ex = :($ex - $bβββ)
end
Expr(:block, :(t = $(esc(x))), :(t2 = 2t), ex)
end
```

Let's try a simple test case: performing Chebyshev interpolation of $\exp(x)$:

```
julia> c = chebcoef(exp)
julia> x = linspace(-1,1,100)
julia> maximum(abs.(Float64[evalcheb(y,c) for y in x] - exp.(x))) # the maximum error on [-1,1]
# check that the evalcheb macro works
julia> evalcheb(0.1234, c[1:4]) - @evalcheb(0.1234, c[1],c[2],c[3],c[4])
```

```
# extract parameters from PowerLawScaling type
julia> pqsPowerLawScaling{p,q,s}(::PowerLawScaling{p,q,s}) = (p,q,s)
```

Extract parameters from PowerLawScaling type

```
julia> pqsPowerLawScaling{p,q,s}(::PowerLawScaling{p,q,s}) = (p,q,s)
julia> @generated function (::FirstIntegral{P,n}, X::Real) where {P<:PowerLawScaling,n}
# compute the Chebyshev coefficients (of the rescaled π¦β as described above)
K = P()
p,q,s = pqsPowerLawScaling(K)
π¦β = X -> quadgk(w -> w^n * K(w*X), 0,1, abstol=1e-12, reltol=1e-10)[1]
Lβ = p < 0 ? X -> π¦β(X) / (s^p + X^p) : π¦β # scale out X βͺ s singularity
q > 0 && throw(DomainError()) # don't know how to deal with growing kernels
qinv = 1/q
c = chebcoef(ΞΎ -> Lβ((1-ΞΎ)^qinv - 2^qinv), 1e-9)
# return an expression that inlines the evaluation of π¦β via C(ΞΎ)
quote
X <= 0 && throw(DomainError())
ΞΎ = 1 - (X + $(2^qinv))^$q
C = @evalcheb ΞΎ $(c...)
return $p < 0 ? C * (X^$p + $(s^p)) : C
end
end
```

A simple example where the result is known analytically:

```
julia> immutable DumbPowerLaw{p,s} <: PowerLawScaling{p,p,s}; end # rα΅ power law
julia> (::FirstIntegral{DumbPowerLaw{p,s}})(r) where {p,s} = r^p
julia> F = FirstIntegral{DumbPowerLaw{-1,1.0},3}()
julia> F(3.7)
julia> @code_llvm F(3.7)
```

```
julia> #Pkg.add("PyPlot")
julia> using PyPlot
julia> x = [0.01:.0125:1.0;];
julia> plot(x, map(FirstIntegral{DumbPowerLaw{-1,1.}, 3}(),x))
```

As an added bonus, I should probably include the other integrals included in section 5.7. Here they are:

Notice that a Helmholtz kernel is also a power law kernel:

```
import GSL
exprel(n, x) = GSL. sf_exprel_n(n, x)
type Helmholtz{k} <: PowerLaw{-1, -1} end # exp(ikr) / 4 pi r
function(::FirstIntegral{Helmholtz{k}, n})(x) where {k,n}
ikx = im * k * x
return exp(ikx) * exprel(n, -ikx)/(4*pi*x)
end
# magnetic field integral equation
abstract type MFIE{k} <: PowerLaw{-3, -3} end # (ikr - 1) * exp(ikr) / 4pi r^3
function (::FirstIntegral{MFIE{k}, n})(x) where {k,n}
ikx = im * k * x
return exp(ikx) * (im * k * exprel(n - 1, -ikx)/((n - 1)*x) - exprel(n - 2, -ikx)/((n-2)*(x^2)) / (4*pi*x)
end
```

© Miguel Raz GuzmΓ‘n Macedo. Last modified: 1 Jan 2000. Website built with Franklin.jl and the Julia programming language.