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?

  • $\begingroup$ 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? $\endgroup$
    – Bort
    Jun 1, 2016 at 8:32

1 Answer 1


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.

These methods are best suited to computing only a small section of the spectrum (for instance, the largest or smallest eigenvalues and their corresponding eigenvectors), and nowadays they are the standard method for solving large and sparse eigenproblems numerically.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.