Evaluation boundary condition hangs for Fourier space solve

[johnbcoughlin]

I am unable to use the Evaluation() operator to set a boundary (gauge) condition on a second-order differential equation in Fourier space.

> ρ = Fun(Fourier(), [0, 0, 1, 2, 3])
> Dx = Derivative(Fourier())

Solving the underspecified Poisson equation works fine:

> Dx^2 \ ρ
Fun(Fourier(【0.0,6.283185307179586❫),[NaN, -0.0, -1.0, -0.5, -0.75])

giving a result with a NaN exactly where you'd expect.

However, trying to set a gauge condition causes the solver to hang:

> [Evaluation(0); Dx^2] \ [0; ρ]
--- uses 100% CPU, appears to be looping ---

Please let me know if I'm maybe using the Evaluation operator incorrectly. The same sort of syntax works for me to solve the examples given in the documentation, and to solve an ODE in terms of Chebyshev functions.

[dlfivefifty]

Just specify a tolerance:

julia> \([Evaluation(0); Dx^2], [0; ρ]; tolerance=100eps())
Fun(Fourier(【0.0,6.283185307179586❫),[1.75, -0.0, -1.0, -0.5, -0.7499999999999999, -2.1912802922805888e-32, -4.3825605845611775e-32, -1.2325951644078312e-32, -8.32667268468867e-17])

If you want to know why, it's because Dx^2 is diagonal:

julia> Dx^2
DerivativeWrapper : Fourier(【0.0,6.283185307179586❫) → Fourier(【0.0,6.283185307179586❫)
 0.0    ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅   -1.0    ⋅     ⋅     ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅   -1.0    ⋅     ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅   -4.0    ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅   -4.0    ⋅     ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅   -9.0    ⋅      ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   -9.0     ⋅      ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅   -16.0     ⋅      ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅   -16.0     ⋅   ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅      ⋅   -25.0  ⋅
  ⋅     ⋅     ⋅     ⋅     ⋅     ⋅     ⋅      ⋅      ⋅      ⋅   ⋱

So adding an extra row destroys the diagonal dominance, which breaks the adaptivity. A numerically better solution is to drop the zero entry:

julia> Dx2 = (Dx^2)[2:end,2:end];

julia> Fun(Dx2 \ ρ, Fourier())
Fun(Fourier(【0.0,6.283185307179586❫),[0.0, -0.0, -1.0, -0.5, -0.75])

[johnbcoughlin]

Thank you, that is helpful for knowing how to use the library. Is it feasible to detect this before attempting the solve routine, or to bail out after some large number of iterations? It is not the nicest experience for library consumers. I am happy to attempt a fix if you can point me to the right location.

[dlfivefifty]

There is a max length of a million. The issue is that building the operator is too slow, for no good reason... this is being redesigned

You can specify a maxlength like:

julia> \([Evaluation(0); Dx^2], [0; ρ]; maxlength=10)
┌ Warning: Maximum length 10 reached.
└ @ ApproxFunBase ~/.julia/packages/ApproxFunBase/Gra7W/src/Caching/matrix.jl:124
Fun(Fourier(【0.0,6.283185307179586❫),[1.75, -0.0, -1.0, -0.5, -0.7499999999999999, -2.1912802922805888e-32, -4.3825605845611775e-32, -1.2325951644078312e-32, -1.2325951644078312e-32, 7.888609052210117e-33, -5.3290705182007475e-17])