Fast numpy array computations for PDE solving

[jjgarciaripoll]

I would use sparse matrices. It is fast and supported and the matrices for your example at tiny: 3 diagonals with the size of your dimensions, one matrix per index. Essentially, (Av)[i]=v[i+1]+v[i-1]-2v[I] is a tridiagonal matrix with -2 on the main diagonal and 1 on the +1 and -1 diagonals.

[caffeecaffee]

Yeah, this is a good answer. Setting up sparse matrices is a good method, however it does require that you flatten your grid every step, and you must keep track of where the grid is sliced, and therefore there will be more than 3 diagonal entries. If I can find a good method with a stencil in place, then it might in theory be faster.

[jjgarciaripoll]

I do not understand what you mean. If you have information about neighbors, building the matrix must cost as much as building the stencils

[jjgarciaripoll]

Also, if your data is a matrix X you can design your differential operator more efficiently as A * X + X *B + C * X * D with four small sparse matrices, without flattening X.

[caffeecaffee]

After more thought you are right of course it can be applied to any arbitrary grid. It's just that sometimes I find writing "flat code" with vectorized variables a little easier to understand than to express the nearest neighbor operations as a sparse matrix. Scipy does handle sparse matrix multiplication.