Integration on an annulus

In this example, we explore integration of the function:

\[ f(x,y) = \frac{x^3}{x^2+y^2-\frac{1}{4}},\]

over the annulus defined by $\{(r,\theta) : \rho < r < 1, 0 < \theta < 2\pi\}$ with parameter $\rho = \frac{2}{3}$. We will calculate the integral:

\[ \int_0^{2\pi}\int_{\frac{2}{3}}^1 f(r\cos\theta,r\sin\theta)^2r{\rm\,d}r{\rm\,d}\theta,\]

by analyzing the function in an annulus polynomial series. We analyze the function on an $N\times M$ tensor product grid defined by:

\[\begin{aligned} r_n & = \sqrt{\cos^2\left[(n+\tfrac{1}{2})\pi/2N\right] + \rho^2 \sin^2\left[(n+\tfrac{1}{2})\pi/2N\right]},\quad{\rm for}\quad 0\le n < N,\quad{\rm and}\\ \theta_m & = 2\pi m/M,\quad{\rm for}\quad 0\le m < M; \end{aligned}\]

we convert the function samples to Chebyshev×Fourier coefficients using plan_annulus_analysis; and finally, we transform the Chebyshev×Fourier coefficients to annulus polynomial coefficients using plan_ann2cxf.

For the storage pattern of the arrays, please consult the documentation.

using FastTransforms, LinearAlgebra, Plots
const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated")
!isdir(GENFIGS) && mkdir(GENFIGS)
plotlyjs()
Plots.PlotlyJSBackend()

Our function $f$ on the annulus:

f = (x,y) -> x^3/(x^2+y^2-1/4)
#1 (generic function with 1 method)

The annulus polynomial degree:

N = 8
M = 4N-3
29

The annulus inner radius:

ρ = 2/3
0.6666666666666666

The radial grid:

r = [begin t = (N-n-0.5)/(2N); ct = sinpi(t); st = cospi(t); sqrt(ct^2+ρ^2*st^2) end for n in 0:N-1]
8-element Vector{Float64}:
 0.9973277184004194
 0.9763124517373388
 0.9362410410518701
 0.8811435628419545
 0.8173313074308551
 0.7535895152498838
 0.7008983100472356
 0.6706577864713554

The angular grid (mod $\pi$):

θ = (0:M-1)*2/M
0.0:0.06896551724137931:1.9310344827586206

On the mapped tensor product grid, our function samples are:

F = [f(r*cospi(θ), r*sinpi(θ)) for r in r, θ in θ]
8×29 Matrix{Float64}:
 1.33215  1.24089  0.995869  0.672118  0.361447  0.136908  0.0255072  0.000211389  -0.00564085  -0.0675532  -0.235438  -0.509748  -0.838068  -1.13371  -1.30886  -1.30886  -1.13371  -0.838068  -0.509748  -0.235438  -0.0675532  -0.00564085  0.000211389  0.0255072  0.136908  0.361447  0.672118  0.995869  1.24089
 1.32342  1.23275  0.989337  0.66771   0.359076  0.13601   0.0253399  0.000210003  -0.00560385  -0.0671101  -0.233894  -0.506405  -0.832572  -1.12628  -1.30028  -1.30028  -1.12628  -0.832572  -0.506405  -0.233894  -0.0671101  -0.00560385  0.000210003  0.0253399  0.13601   0.359076  0.66771   0.989337  1.23275
 1.30981  1.22008  0.979168  0.660847  0.355385  0.134612  0.0250795  0.000207844  -0.00554625  -0.0664203  -0.23149   -0.5012    -0.824014  -1.1147   -1.28691  -1.28691  -1.1147   -0.824014  -0.5012    -0.23149   -0.0664203  -0.00554625  0.000207844  0.0250795  0.134612  0.355385  0.660847  0.979168  1.22008
 1.29961  1.21057  0.97154   0.655698  0.352617  0.133563  0.0248841  0.000206225  -0.00550305  -0.0659028  -0.229687  -0.497295  -0.817594  -1.10601  -1.27689  -1.27689  -1.10601  -0.817594  -0.497295  -0.229687  -0.0659028  -0.00550305  0.000206225  0.0248841  0.133563  0.352617  0.655698  0.97154   1.21057
 1.30613  1.21665  0.976415  0.658989  0.354386  0.134233  0.025009   0.00020726   -0.00553066  -0.0662336  -0.230839  -0.499791  -0.821697  -1.11156  -1.28329  -1.28329  -1.11156  -0.821697  -0.499791  -0.230839  -0.0662336  -0.00553066  0.00020726   0.025009   0.134233  0.354386  0.658989  0.976415  1.21665
 1.34623  1.25399  1.00639   0.679218  0.365265  0.138354  0.0257767  0.000213622  -0.00570044  -0.0682668  -0.237925  -0.515133  -0.846922  -1.14569  -1.32269  -1.32269  -1.14569  -0.846922  -0.515133  -0.237925  -0.0682668  -0.00570044  0.000213622  0.0257767  0.138354  0.365265  0.679218  1.00639   1.25399
 1.42719  1.32941  1.06692   0.720069  0.387234  0.146675  0.027327   0.00022647   -0.00604329  -0.0723726  -0.252235  -0.546115  -0.897858  -1.21459  -1.40224  -1.40224  -1.21459  -0.897858  -0.546115  -0.252235  -0.0723726  -0.00604329  0.00022647   0.027327   0.146675  0.387234  0.720069  1.06692   1.32941
 1.5099   1.40645  1.12874   0.761795  0.409673  0.155175  0.0289105  0.000239594  -0.00639348  -0.0765664  -0.266852  -0.577762  -0.949887  -1.28498  -1.4835   -1.4835   -1.28498  -0.949887  -0.577762  -0.266852  -0.0765664  -0.00639348  0.000239594  0.0289105  0.155175  0.409673  0.761795  1.12874   1.40645

We superpose a surface plot of $f$ on top of the grid:

X = [r*cospi(θ) for r in r, θ in θ]
Y = [r*sinpi(θ) for r in r, θ in θ]
scatter3d(vec(X), vec(Y), vec(0F); markersize=0.75, markercolor=:red)
surface!(X, Y, F; legend=false, xlabel="x", ylabel="y", zlabel="f")
savefig(joinpath(GENFIGS, "annulus.html"))
"/juliateam/.julia/packages/FastTransforms/0VwwO/docs/src/generated/annulus.html"

We precompute an Annulus–Chebyshev×Fourier plan:

α, β, γ = 0, 0, 0
P = plan_ann2cxf(F, α, β, γ, ρ)
FastTransforms Annulus--Chebyshev×Fourier plan for 8×29-element array of Float64

And an FFTW Chebyshev×Fourier analysis plan on the annulus:

PA = plan_annulus_analysis(F, ρ)
FastTransforms plan for FFTW Chebyshev×Fourier analysis on the annulus for 8×29-element array of Float64

Its annulus coefficients are:

U = P\(PA*F)
8×29 Matrix{Float64}:
 -2.45849e-17   3.1221e-17    0.926725      2.13898e-17  -7.722e-18     6.63286e-17   0.29418      -8.04495e-18   2.67464e-17   2.18818e-17  -6.62939e-18   1.7977e-17   -4.74235e-18   2.27738e-18  -3.14548e-17   5.7696e-17    2.61652e-17   2.19612e-17   7.8769e-18    3.2266e-17    7.7026e-17   -1.27667e-17   5.90664e-17  -1.50384e-17  -3.21534e-17  -4.94481e-18   3.24342e-17  -3.0763e-17   -5.02811e-17
 -9.36993e-18  -5.95366e-18  -0.124037     -1.22929e-17   3.70137e-17  -2.70714e-17  -0.0976704    -8.08942e-18  -2.32787e-17  -1.69012e-17   2.00575e-17  -1.47269e-17  -1.34007e-17  -7.29501e-18   3.03374e-17  -3.68364e-17  -1.91986e-17  -1.70503e-17  -8.74329e-18  -1.78914e-17  -6.71269e-17   1.39282e-17  -5.51026e-17   1.76125e-17   1.75724e-17   1.23183e-17  -3.17445e-17   2.56533e-17   5.34817e-17
 -2.1258e-17   -1.71836e-18   0.042116     -3.07069e-18  -1.44328e-17   6.118e-18     0.0336531    -4.96328e-19   1.41794e-17   7.80795e-18  -9.39246e-18   7.39772e-18   3.20181e-18   1.87799e-18  -1.35943e-17   2.18868e-17   6.20453e-18   6.34172e-18  -2.37815e-19   1.08571e-17   2.30856e-17  -1.57391e-17   3.67156e-17  -1.92223e-17  -3.74272e-17  -9.23703e-18   1.82113e-17  -2.24539e-17  -5.23781e-17
 -5.78545e-18   4.26889e-18  -0.0139459     5.99855e-18  -1.97145e-18   1.06084e-17  -0.011224      1.07487e-17   2.25101e-17   9.78784e-18  -1.706e-18     8.31177e-18   8.92314e-18   9.3792e-18    1.05443e-17   3.29865e-18  -3.8951e-18    3.61053e-18   5.94666e-18  -4.5099e-19   -1.87384e-17   4.16053e-18  -1.97598e-17   2.12786e-18   7.40023e-18  -6.73183e-19  -1.61549e-17   6.4625e-18    4.61626e-18
 -2.40148e-17   5.3755e-18    0.00458845    3.58419e-18   1.86344e-17   6.63249e-18   0.00369818   -1.96889e-18  -1.04804e-17  -2.72812e-18   1.0092e-17   -5.36427e-18   6.14311e-18  -7.10032e-18  -3.07984e-18  -2.05996e-18   4.32616e-19  -3.88475e-18  -1.10073e-18  -9.90196e-19   1.25448e-17  -4.46876e-18   1.48269e-17  -3.54572e-18   4.78504e-18  -2.09469e-18   1.21511e-17  -3.97586e-18   5.39528e-18
  2.61466e-17   1.24044e-19  -0.00150441   -1.05751e-18  -9.38596e-18  -5.84846e-18  -0.00121128   -2.20622e-18  -8.79109e-18  -5.36301e-18   2.65743e-18  -4.71053e-18   7.8515e-18   -2.78458e-18   8.25058e-18  -6.29463e-18  -5.34413e-18  -3.77451e-19   1.52976e-18   5.26751e-20  -2.3092e-19    6.15042e-18   1.45697e-18   7.07569e-18   9.27665e-18   4.31379e-18   6.65377e-18   5.57142e-18   4.29955e-18
  6.05256e-17  -5.88287e-18   0.000517281  -3.58556e-18   1.29799e-17  -8.61463e-18   0.000424358  -2.99387e-18   3.91183e-18  -3.46929e-18  -9.8607e-18    2.57698e-19  -3.43595e-18   8.66221e-19  -1.85307e-17  -1.59528e-18   5.24491e-18   2.68922e-18  -1.04194e-17  -2.44224e-18  -2.51142e-18   1.0758e-18   -1.09894e-17   1.03318e-18   1.59905e-18  -5.57554e-19   1.2693e-18    1.49178e-18   5.50016e-18
 -1.16089e-17   1.73822e-18  -0.000162126   3.61748e-18   7.44168e-18   3.44124e-18  -5.40419e-5    3.85381e-18   7.08965e-18   2.94047e-18   6.55323e-18   2.13461e-18  -9.07416e-18   7.16234e-19  -3.88089e-19  -7.35244e-19  -1.12942e-17  -1.94874e-18   1.19628e-17  -2.90149e-18  -1.04048e-17  -3.459e-18     1.06577e-17  -3.54894e-18  -1.70683e-17  -2.75522e-18   4.99216e-18  -1.13437e-18  -7.73268e-18

The annulus coefficients are useful for integration. The integral of $[f(x,y)]^2$ over the annulus is approximately the square of the 2-norm of the coefficients:

norm(U)^2, 5π/8*(1675/4536+9*log(3)/32-3*log(7)/32)
(0.9735516844404255, 0.973547572736036)

This page was generated using Literate.jl.