Poisson Solvers

Solve

\[- \Delta \phi(x) = f(x)\]

using Plots
using PoissonSolvers

sol(x) = sin(π*x) * (1-x^2)
der(x) = π * cos(π*x) * (1-x^2) - 2 * x * sin(π*x)
rhs(x) = 4π * x * cos(π*x) + 2 * sin(π*x) + π^2 * (1-x^2) * sin(π*x)

domain = (-1., +1.)

x = LinRange(domain[begin], domain[end], 100)
100-element LinRange{Float64, Int64}:
 -1.0, -0.979798, -0.959596, -0.939394, …, 0.939394, 0.959596, 0.979798, 1.0

FFT Solver

b = FFTWBasis(domain, 256)
p = Potential(b, rhs)
Potential{FFTWBasis{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{Float64}, PoissonSolvers.FFTWSolution{Vector{Float64}, FFTWBasis{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, PoissonSolverFFT{Float64, FFTWBasis{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}(FFTWBasis{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}((-1.0, 1.0), [..., -1.0, -0.9921875, -0.984375, -0.9765625, -0.96875, -0.9609375, -0.953125, -0.9453125, -0.9375, -0.9296875  …  0.9296875, 0.9375, 0.9453125, 0.953125, 0.9609375, 0.96875, 0.9765625, 0.984375, 0.9921875, 1.0, ...], 0.0078125), [0.016284748572368737, 0.015539437902630574, 0.014015151427750776, 0.011756235249516911, 0.008759460926567382, 0.005044954857989926, 0.0006205622830512958, -0.0044964107492057315, -0.010293982587165157, -0.016754817222198015  …  0.04921868670067725, 0.04281989249444884, 0.037074245745341405, 0.03199983751570473, 0.027609427762632865, 0.023921124366353796, 0.02094359572330967, 0.018697814708024174, 0.017181416293134577, 0.016439641820030215], PoissonSolvers.FFTWSolution{Vector{Float64}, FFTWBasis{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}(FFTWBasis{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}((-1.0, 1.0), [..., -1.0, -0.9921875, -0.984375, -0.9765625, -0.96875, -0.9609375, -0.953125, -0.9453125, -0.9375, -0.9296875  …  0.9296875, 0.9375, 0.9453125, 0.953125, 0.9609375, 0.96875, 0.9765625, 0.984375, 0.9921875, 1.0, ...], 0.0078125), [0.016284748572368737, 0.015539437902630574, 0.014015151427750776, 0.011756235249516911, 0.008759460926567382, 0.005044954857989926, 0.0006205622830512958, -0.0044964107492057315, -0.010293982587165157, -0.016754817222198015  …  0.04921868670067725, 0.04281989249444884, 0.037074245745341405, 0.03199983751570473, 0.027609427762632865, 0.023921124366353796, 0.02094359572330967, 0.018697814708024174, 0.017181416293134577, 0.016439641820030215]), PoissonSolverFFT{Float64, FFTWBasis{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}(FFTWBasis{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}((-1.0, 1.0), [..., -1.0, -0.9921875, -0.984375, -0.9765625, -0.96875, -0.9609375, -0.953125, -0.9453125, -0.9375, -0.9296875  …  0.9296875, 0.9375, 0.9453125, 0.953125, 0.9609375, 0.96875, 0.9765625, 0.984375, 0.9921875, 1.0, ...], 0.0078125)), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
plot(xlabel = "x", ylabel = "ϕ(x)")
plot!(x, p.(x); xlims = domain, label = "Solution")
plot!(x, sol.(x); xlims = domain, label = "Reference")
Example block output
plot(xlabel = "x", ylabel = "ϕ'(x)")
plot!(x, p.(x, Derivative(1)); xlims = domain, label = "Derivative")
plot!(x, der.(x); xlims = domain, label = "Reference")
Example block output

B-Spline Solver

b = PeriodicBasisBSplineKit(domain, 3, 32)
p = Potential(b, rhs)
Potential{BSplineKit.BSplines.PeriodicBSplineBasis{3, Float64, BSplineKit.BSplines.PeriodicKnots{Float64, 3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, Vector{Float64}, BSplineKit.Splines.Spline{Float64, BSplineKit.BSplines.PeriodicBSplineBasis{3, Float64, BSplineKit.BSplines.PeriodicKnots{Float64, 3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, BSplineKit.Splines.PeriodicVector{Float64, Vector{Float64}}}, PoissonSolverBSplineKit{Float64, ComplexF64, BSplineKit.BSplines.PeriodicBSplineBasis{3, Float64, BSplineKit.BSplines.PeriodicKnots{Float64, 3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}}(32-element PeriodicBSplineBasis of order 3, domain [-1.0, 1.0), period 2.0
 knots: [..., -1.0625, -1.0, -0.9375, -0.875, -0.8125, -0.75, -0.6875, -0.625, -0.5625, -0.5  …  0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1.0, 1.0625, ...], [0.01363183816829476, -0.041278025565521355, -0.12503498107439626, -0.23903728465449753, -0.3675433153209978, -0.49816096746892635, -0.6176180884241456, -0.7141028308697789, -0.7777151169942914, -0.8012307571265946  …  0.8012307571265948, 0.7777151169942913, 0.7141028308697791, 0.6176180884241455, 0.4981609674689263, 0.3675433153209976, 0.2390372846544975, 0.1250349810743961, 0.04127802556552118, -0.013631838168294948], 32-element Spline{Float64}:
 basis: 32-element PeriodicBSplineBasis of order 3, domain [-1.0, 1.0), period 2.0
 order: 3
 knots: [..., -1.0625, -1.0, -0.9375, -0.875, -0.8125, -0.75, -0.6875, -0.625, -0.5625, -0.5  …  0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1.0, 1.0625, ...]
 coefficients: [..., 0.0136318, -0.041278, -0.125035, -0.239037, -0.367543, -0.498161, -0.617618, -0.714103, -0.777715, -0.801231  …  0.801231, 0.777715, 0.714103, 0.617618, 0.498161, 0.367543, 0.239037, 0.125035, 0.041278, -0.0136318, ...], PoissonSolverBSplineKit{Float64, ComplexF64, BSplineKit.BSplines.PeriodicBSplineBasis{3, Float64, BSplineKit.BSplines.PeriodicKnots{Float64, 3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}}(32-element PeriodicBSplineBasis of order 3, domain [-1.0, 1.0), period 2.0
 knots: [..., -1.0625, -1.0, -0.9375, -0.875, -0.8125, -0.75, -0.6875, -0.625, -0.5625, -0.5  …  0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1.0, 1.0625, ...], [0.034374999999999996 0.013541666666666664 … 0.0005208333333333339 0.013541666666666664; 0.013541666666666664 0.034374999999999996 … 0.0 0.0005208333333333339; … ; 0.0005208333333333339 0.0 … 0.034374999999999996 0.013541666666666664; 0.013541666666666664 0.0005208333333333339 … 0.013541666666666664 0.034374999999999996], [15.999999999999986 -5.3333333333333215 … -2.6666666666666696 -5.3333333333333215; -5.3333333333333215 15.999999999999986 … 0.0 -2.6666666666666696; … ; -2.6666666666666696 0.0 … 15.999999999999986 -5.3333333333333215; -5.3333333333333215 -2.6666666666666696 … -5.3333333333333215 15.999999999999986], ToeplitzMatrices.ToeplitzFactorization{Float64, ToeplitzMatrices.Circulant{Float64, SparseArrays.SparseVector{Float64, Int64}}, ComplexF64, FFTW.cFFTWPlan{ComplexF64, -1, true, 1, Tuple{Int64}}}(ComplexF64[0.06249999999999999 + 0.0im, 0.061900309190620076 + 0.0im, 0.060133306902583325 + 0.0im, 0.057292597241907564 + 0.0im, 0.05352580865713565 + 0.0im, 0.04902306523556725 + 0.0im, 0.044002773396151856 + 0.0im, 0.038696321708237544 + 0.0im, 0.033333333333333326 + 0.0im, 0.028128929265697267 + 0.0im  …  0.02327408747637616 + 0.0im, 0.028128929265697267 + 0.0im, 0.033333333333333326 + 0.0im, 0.038696321708237544 + 0.0im, 0.044002773396151856 + 0.0im, 0.04902306523556725 + 0.0im, 0.05352580865713565 + 0.0im, 0.057292597241907564 + 0.0im, 0.060133306902583325 + 0.0im, 0.061900309190620076 + 0.0im], ComplexF64[5.0e-324 + 3.5e-323im, 6.92099716539175e-310 + 6.92099716539333e-310im, 4.0e-323 + 4.0e-323im, 6.9209971653949e-310 + 6.9209971653965e-310im, 4.4e-323 + 4.4e-323im, 6.92099716539808e-310 + 6.92099716539966e-310im, 5.0e-323 + 7.0e-323im, 6.92099716540124e-310 + 6.9209971654028e-310im, 7.4e-323 + 9.0e-323im, 6.9209971654044e-310 + 6.920997165406e-310im  …  1.9e-322 + 1.93e-322im, 6.92099716542653e-310 + 6.9209971654281e-310im, 2.0e-322 + 2.17e-322im, 6.9209971654297e-310 + 6.9209971654313e-310im, 2.2e-322 + 2.2e-322im, 6.92099716543286e-310 + 6.92099716543444e-310im, 2.27e-322 + 2.27e-322im, 6.920997165436e-310 + 6.9209971654376e-310im, 6.92104469222944e-310 + 2.5e-323im, 3.0e-323 + 3.5e-323im], FFTW in-place forward plan for 32-element array of ComplexF64
(dft-direct-32 "n2fv_32_avx2_128")), ToeplitzMatrices.ToeplitzFactorization{Float64, ToeplitzMatrices.Circulant{Float64, SparseArrays.SparseVector{Float64, Int64}}, ComplexF64, FFTW.cFFTWPlan{ComplexF64, -1, true, 1, Tuple{Int64}}}(ComplexF64[1.0000000000000036 + 0.0im, 0.6109328356386836 + 0.0im, 2.374048820218025 + 0.0im, 5.0900124961590425 + 0.0im, 8.457527667343495 + 0.0im, 12.114895820404726 + 0.0im, 15.689279554433963 + 0.0im, 18.846394071888156 + 0.0im, 21.333333333333325 + 0.0im, 23.00832094156555 + 0.0im  …  23.853192778222525 + 0.0im, 23.00832094156555 + 0.0im, 21.333333333333325 + 0.0im, 18.846394071888156 + 0.0im, 15.689279554433963 + 0.0im, 12.114895820404726 + 0.0im, 8.457527667343495 + 0.0im, 5.0900124961590425 + 0.0im, 2.374048820218025 + 0.0im, 0.6109328356386836 + 0.0im], ComplexF64[0.01363183816829476 + 0.0im, -0.041278025565521355 + 0.0im, -0.12503498107439626 + 0.0im, -0.23903728465449753 - 3.469446951953614e-18im, -0.3675433153209978 + 0.0im, -0.49816096746892635 - 3.469446951953614e-18im, -0.6176180884241456 + 0.0im, -0.7141028308697789 + 0.0im, -0.7777151169942914 + 0.0im, -0.8012307571265946 + 0.0im  …  0.8012307571265948 + 0.0im, 0.7777151169942913 + 0.0im, 0.7141028308697791 + 0.0im, 0.6176180884241455 + 0.0im, 0.4981609674689263 + 0.0im, 0.3675433153209976 + 3.469446951953614e-18im, 0.2390372846544975 + 0.0im, 0.1250349810743961 + 3.469446951953614e-18im, 0.04127802556552118 + 0.0im, -0.013631838168294948 + 0.0im], FFTW in-place forward plan for 32-element array of ComplexF64
(dft-direct-32 "n2fv_32_avx2_128")), [0.96875 -0.03125 … -0.03125 -0.03125; -0.03125 0.96875 … -0.03125 -0.03125; … ; -0.03125 -0.03125 … 0.96875 -0.03125; -0.03125 -0.03125 … -0.03125 0.96875], [0.03125 0.03125 … 0.03125 0.03125; 0.03125 0.03125 … 0.03125 0.03125; … ; 0.03125 0.03125 … 0.03125 0.03125; 0.03125 0.03125 … 0.03125 0.03125]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
plot(xlabel = "x", ylabel = "ϕ(x)")
plot!(x, p.(x); xlims = domain, label = "Solution")
plot!(x, sol.(x); xlims = domain, label = "Reference")
Example block output
plot(xlabel = "x", ylabel = "ϕ'(x)")
plot!(x, p.(x, Derivative(1)); xlims = domain, label = "Derivative")
plot!(x, der.(x); xlims = domain, label = "Reference")
Example block output