Solve eigenvalue problem using finite differences without vectorization

I am interested in solving the problem $$-A u = \lambda u$$ on a finite differences grid on a square. In my case, the operator $$A$$ is of the type $$-\Delta + \mu I$$, where $$\mu$$ is there to impose some constraints. The method I used until now was to consider a finite differences grid and get $$4u_{i,j} - u_{i,j+1}-u_{i,j-1}-u_{i+1,j}-u_{i-1,j} + \mu_{i,j} u_{i,j}=\lambda u_{i,j}$$ Then we transform $$u_{i,j}$$ in a vector and write the above problem in matrix form $$Au = \lambda u$$. The problem I have with this is that the size of the matrix $$A$$ is $$n \times n$$, where $$n$$ is the number of points in the grid. Therefore, for a grid with $$N^2$$ points the dimension of the matrix is $$N^2 \times N^2$$. Therefore, the size of the matrices grows very fast with respect to the discretization.

Is it possible to solve this kind of problem without vectorization, i.e. using only matrices of size $$N \times N$$?

My idea is the following: If we denote $$U=(u_{i,j})$$ then the discrete PDE can be written as $$(4I -S_{0,1}-S_{0,-1}-S_{1,0}-S_{-1,0})U+\mu \odot U=\lambda U,$$ where $$S_{i,j}$$ is the shift of the matrix with $$i$$ on horizontal and $$j$$ on vertical directions. The $$\odot$$ is "pointwise" product. Can this lead to a solution?

• Given your discretization scheme, I would assume that you should be able to use large N with a sparse matrix represention and sparse matrix solver. Have you considered those? – Bort Jun 1 '16 at 8:32

For sure it is possible with iterative methods, because you can apply your operator $u\mapsto Au$ cheaply while keeping $u$ in "matrix form". It is simply a double nested for cycle. Check out ARPACK, or more in general Arnoldi-(or Krylov-)type methods: in most implementations, you can use as input a 'black-box' routine that computes $Au$ from $u$, rather than an explicit matrix.