Integral operators (again)

variational formulation
3

So there was certainly a finite error in the symmetry breaking, and it did resolve itself under mesh refinement. Your suggestion of symmetrizing the matrix is a decent one. The issue is that, effectively, your implementation treats the two arguments in the kernel different; one argument is projected onto the FEM space, while the other argument is interpolated onto it at the nodal values.

I came up with an alternative, which is to treat it via interpolation in both arguments, so that

c(x,y) \approx \sum_{i,j} c(x_i,y_j) u_i(x) u_j(y)

which induces a matrix C, with C_{i,j} = c(x_i,y_j). Then, you have the symmetric linear system:

M C M \phi= \lambda M \phi

and this provides satisfactory results. This may also be of use in the Fredholm equation that you were originally interestsed.

Two issues I am still less than thrilled about, but may be unavoidable are as follows:

  1. I am using dense NumPy matrices for C because I’m not sure there’s a better option. My kernel does fall off relatively quickly with |x-y|, so maybe if I introduced a cutoff I could justify using a sparse representation.
  2. Since I am using these dense matrix representations, I end up just using the SciPy eigh eigensolver. Does anyone know if the SLEPc dense solvers would be more efficient, were I to set the problem up with them?