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/CuqEn/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   2.92896e-17   0.926725      2.04896e-17  -7.722e-18     6.1888e-17    0.29418      -1.81602e-17   2.67464e-17   1.20498e-17  -6.62939e-18   1.11857e-17  -1.67703e-17  -6.13512e-19  -3.42122e-17   6.23458e-17   1.8968e-17    3.12444e-17   4.92041e-18   3.83711e-17   7.86308e-17  -1.08653e-17   5.78808e-17  -1.56088e-17  -3.34631e-17  -8.2245e-18    3.76837e-17  -3.30615e-17  -4.36512e-17
 -9.36993e-18  -2.254e-18    -0.124037     -5.91733e-18   3.70137e-17  -1.89099e-17  -0.0976704     1.024e-18    -2.32787e-17  -9.77876e-18   2.00575e-17  -1.04302e-17   7.97476e-19  -5.50321e-18   4.69595e-17  -4.01755e-17  -1.55142e-18  -2.52689e-17   1.02515e-18  -2.67731e-17  -6.90711e-17   5.4336e-18   -5.62647e-17   1.0013e-17    1.53665e-17   7.96412e-18  -4.27121e-17   2.50418e-17   4.09539e-17
 -2.1258e-17    1.12715e-18   0.042116      9.94105e-19  -1.44328e-17   9.15449e-18   0.0336531     3.6262e-19    1.41794e-17   7.39131e-18  -9.39246e-18   6.73552e-18   7.9017e-18    1.31493e-18  -1.74808e-17   2.30498e-17  -4.92523e-18   1.03454e-17  -1.03774e-17   1.67754e-17   1.97933e-17  -8.57e-18      3.57962e-17  -1.14124e-17  -3.35749e-17  -5.45601e-18   2.98569e-17  -2.14064e-17  -3.95395e-17
 -5.78545e-18   4.00391e-18  -0.0139459     3.14722e-18  -1.97145e-18   5.82968e-18  -0.011224      6.26797e-18   2.25101e-17   6.88701e-18  -1.706e-18     7.40166e-18   3.44473e-18   9.55486e-18   5.33772e-18   2.9721e-18   -8.00488e-18   3.40533e-18   4.2453e-18    2.2961e-19   -1.44173e-17   7.73146e-18  -2.06172e-17   5.91455e-18   1.1306e-17    3.25478e-18  -1.21454e-17   8.18752e-18   4.67076e-18
 -2.40148e-17   7.71573e-19   0.00458845   -2.03822e-18   1.86344e-17   3.19547e-18   0.00369818   -1.86152e-18  -1.04804e-17   8.11326e-19   1.0092e-17   -1.77139e-18  -1.80973e-18  -6.8012e-18    1.80129e-18  -4.24712e-18   6.46529e-18  -7.2048e-18    7.11746e-18  -4.72462e-18   1.73345e-17  -5.27319e-18   1.40121e-17  -1.67467e-18  -4.79598e-19   8.00009e-19   6.60738e-18  -3.34661e-18  -2.0813e-18
  2.61466e-17   1.71852e-19  -0.00150441    7.97602e-19  -9.38596e-18  -3.21855e-18  -0.00121128   -4.47164e-19  -8.79109e-18  -4.83245e-18   2.65743e-18  -5.47568e-18  -5.50836e-19  -2.9175e-18    7.42411e-18  -5.70705e-18  -3.31104e-18  -3.85137e-19   2.72442e-18  -5.31797e-19  -2.48579e-18   3.59849e-18   7.3637e-19    3.4972e-18    6.6109e-18    1.53405e-18   6.20867e-18   4.33464e-18   7.06298e-18
  6.05256e-17  -4.20047e-18   0.000517281  -2.62968e-18   1.29799e-17  -9.06153e-18   0.000424358  -4.57265e-18   3.91183e-18  -6.80939e-18  -9.8607e-18   -2.399e-18     2.87152e-18   6.6648e-19   -2.13866e-17   1.69496e-20   4.87425e-18   5.37575e-18  -1.38669e-17   3.57204e-19  -5.42692e-18   2.48838e-18  -1.15939e-17   4.39475e-19   3.73151e-18  -1.63971e-18   2.61275e-18   1.39577e-18   8.50526e-18
 -1.16089e-17   6.39707e-19  -0.000162126   3.55753e-18   7.44168e-18   4.28763e-18  -5.40419e-5    5.16629e-18   7.08965e-18   5.4486e-18    6.55323e-18   4.64571e-18   2.63303e-18   6.04678e-19   6.15346e-19  -2.13942e-18  -1.32656e-17  -3.57457e-18   1.16079e-17  -5.67995e-18  -1.01101e-17  -4.98049e-18   1.02985e-17  -2.89642e-18  -1.8487e-17   -2.30036e-18   2.95997e-18  -1.26575e-18  -1.033e-17

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.