Julia package for computing multivariate singular integrals
The function newtoniansquare([x,y], p)
computes a matrix of Newtonian potentials of Legendre polynomials on the unit square
where
where
Here we show an example of how to use the package for computing the Newtonian potential for
julia> using MultivariateSingularIntegrals, ClassicalOrthogonalPolynomials, LinearAlgebra, QuadGK
julia> f = (x,y) -> cos(x * exp(y));
julia> 𝐱 = [0.1, 1.1]; # point near the unit square
julia> @time qgk_approx = quadgk(s -> quadgk(t -> f(s, t) * log(norm(𝐱 - [s,t])), -1, 1)[1], -1, 1)[1]; # compute the integral using quadgk
0.388834 seconds (804.74 k allocations: 35.045 MiB, 3.99% gc time, 98.86% compilation time)
julia> P = Legendre(); # Legendre polynomials
julia> p = 40; # degree of the Legendre polynomials
julia> x = ClassicalOrthogonalPolynomials.grid(P, p); # grid of points
julia> F = f.(x, x') # evaluate f on the grid;
julia> @time C = plan_transform(P, (p, p)) * F; # 2D Legendre coefficients
0.002116 seconds (2.87 k allocations: 3.594 MiB)
julia> @time N = Float64.(newtoniansquare(big.(𝐱), p)); # compute the Newtonian potential of Legendre polynomials, using BigFloat to avoid numerical issues
0.003480 seconds (80.94 k allocations: 4.119 MiB)
julia> @test dot(N, C) ≈ qgk_approx
Test Passed
The package continues to work inside the unit square:
julia> 𝐱 = [0.1, 0.2];
julia> @time N = Float64.(newtoniansquare(big.(𝐱), p));
0.005985 seconds (165.52 k allocations: 8.324 MiB)
julia> dot(N, C)
-1.2555132824835833
The package also supports complex logarithmic integrals with kernel logkernelsquare
and Stieltjes integrals with kernel stieltjessquare
. The real and imaginary parts of Stieltjes integrals containing the gradient of the Newtonian potential.