Support zero Dirichlet conditions

[dlfivefifty]

No description provided.

[dlfivefifty]

@ioannisPApapadopoulos @KarsKnook Any ideas what to about the ADI shifts become NAN when p = 160?

julia> J, a, b, c, d, tol
(78, 9.10002282376538e-10, 0.4052847345693519, -0.4052847345693519, -9.10002282376538e-10, 1.0e-15)

julia> ADI_shifts(J, a, b, c, d, tol)
([NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN  …  NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN], [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN  …  NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN])

It looks like the culprit is the elliptic integral but not even using big works:

julia> α
4.453667231591123e8

julia>     K = Elliptic.K(1-1/big(α)^2)
∞

[codecov]

Codecov Report

Patch coverage: 96.91% and project coverage change: -0.39% ⚠️

Comparison is base (e2dcaec) 98.70% compared to head (daefcd0) 98.32%.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #37      +/-   ##
==========================================
- Coverage   98.70%   98.32%   -0.39%     
==========================================
  Files           4        5       +1     
  Lines         386      537     +151     
==========================================
+ Hits          381      528     +147     
- Misses          5        9       +4     

Files Changed	Coverage Δ
src/PiecewiseOrthogonalPolynomials.jl	100.00% <ø> (ø)
src/dirichlet.jl	96.29% <96.29%> (ø)
src/arrowhead.jl	98.63% <97.16%> (-0.35%)	⬇️
src/continuouspolynomial.jl	97.60% <100.00%> (ø)

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[dlfivefifty]

Hmm it just converts a BigFloat back to a Float64:

https://github.com/nolta/Elliptic.jl/blob/7ed7dfbb0df9398920562b27df7d3cf782b0fa76/src/Elliptic.jl#L100

[TSGut]

Have you considered using hypergeometric functions instead? Hypergeometric functions are great. The formula I think is

Small example for sanity check:

julia> α = 4.453667231591123
4.453667231591123

julia> K = Elliptic.K(1-1/α^2)
2.9043516240612575

julia> K2f1 = convert(eltype(α),π)/2*HypergeometricFunctions._₂F₁(one(eltype(α))/2,one(eltype(α))/2,1,1-1/α^2)
2.904351624061259

julia> @test K ≈ K2f1
Test Passed

and now your example:

julia> α = 4.453667231591123e8
4.453667231591123e8

julia> K = Elliptic.K(1-1/big(α)^2)
Inf

julia> K2f1 = convert(eltype(α),π)/2*HypergeometricFunctions._₂F₁(one(eltype(α))/2,one(eltype(α))/2,1,1-1/big(α)^2)
21.3007029588556941903173410911857761893553940237856013580751239264221275085114

It does need BigFloat for your big example and I have not tested whether the answer is actually accurate at all.

[MikaelSlevinsky]

Have u tried the elliptic functions sub module in FastTransforms?

[MikaelSlevinsky]

There's also code for ADI shifts starting at line 2935 in FastTransforms C library banded_source.c

[MikaelSlevinsky]

I guess I should add BigFloat elliptic support

[MikaelSlevinsky]

Your ADI_shifts needs a case for argument unity that uses an asymptotic expansion https://github.com/MikaelSlevinsky/FastTransforms/blob/master/src/banded_source.c#L2937-L2954

[dlfivefifty]

@TSGut Btw its cleaner to just write one(α) instead of one(eltype(α))

[dlfivefifty]

OK it now "works" except the result is inaccurate with p = 160 (and fine with p = 120)! I suspect I need to replace the dn calls with BigFloat-Hypergeometric as well. @TSGut Know of a good formula?

@tomtrogdon @cade-b Maybe you know a good way to compute dn in Julia?

[dlfivefifty]

(At least I can do a bait-and-switch in my talk and use this just to show the timings....)

[KarsKnook]

For $\alpha > 10^7$ there is a specific way to compute the ADI shifts which I did not implement. See appendix A of Dan and Alex's paper

[dlfivefifty]

Thanks @KarsKnook !! I implemented it and it now works for p = 160 (that is 639^2 DOFs, taking 3.8s).

[tomtrogdon]

I've just used Elliptic.jl. But the methods https://dlmf.nist.gov/22.20 on DLMF seem reasonable to implement.

[MikaelSlevinsky]

Except alpha > 1e7 is specific to Float64. In Float32, that would fail for alpha > 1/sqrt(eps(Float32)). So what you really want is to check that the argument of the complete elliptic integral of the first kind isn't equal to one in the current precision.

[dlfivefifty]

I have no intention of solving Poisson on a square with anything other than Float64....