# Matrix to generate finite difference

I'm implementing a finite difference scheme for a 2D PDE problem. I wish to avoid using a loop to generate the finite differences. For instance to generate a 2nd order central difference of u(x,y)_xx, I can multiply u(x,y) by the following:

Is there a nice matrix representation for u_xy = (u_{i+1,j+1} + u_{i-1,j-1} - u_{i-1,j+1} - u_{i+1,j-1})/(4dxdy)? It's a harder problem to code as it's in 2D - I'd like to multiply some matrix by u(x,y) to avoid looping. Many thanks!

## Answers

If your points are stored in a N-by-N matrix then, as you said, left multiplying by your finite difference matrix gives an approximation to the second derivative with respect to u_{xx}. Right-multiplying by the *transpose* of the finite difference matrix is equivalent to an approximation u_{yy}. You can get an approximation to the mixed derivative u_{xy} by left-multiplying **and** right-multiplying by e.g. a central difference matrix

delta_2x = 0 1 0 0 0 -1 0 1 0 0 0 -1 0 1 0 0 0 -1 0 1 0 0 0 -1 0

(then divide by the factor 4*Dx*Dy), so something like

U_xy = 1/(4*Dx*Dy) * delta_2x * U_matrix * delta_2x';

If you cast a matrix as a N^2 vector

U_vec = U_matrix(:);

then these operators can be expressed using a Kronecker product, implemented in MATLAB as kron: We have

A*X*B = kron(B',A)*X(:);

so for your finite difference matrices

U_xy_vec = 1/(4*Dx*Dy)*(kron(delta_2x,delta_2x)*U_vec);

If instead you have an N-by-M matrix U_mat, then left matrix multiplication is equivalent to kron(eye(M),delta_2x_N) and right multiplication to kron(delta_2y_M,eye(N)), where delta_2y_M (delta_2x_N) is the M-by-M (N-by-N) central difference matrix, so the operation is

U_xy_vec = 1/(4*Dx*Dy) * kron(delta_2y_M,delta_2y_N)*U_vec;

Here is an MATLAB code example:

N = 20; M = 30; Dx = 1/N; Dy = 1/M; [Y,X] = meshgrid((1:(M))./(M+1),(1:(N))/(N+1)); % Example solution and mixed derivative (chosen for 0 BCs) U_mat = sin(2*pi*X).*(sin(2*pi*Y.^2)); U_xy = 8*pi^2*Y.*cos(2*pi*X).*cos(2*pi*Y.^2); % Centred finite difference matrices delta_x_N = 1/(2*Dx)*(diag(ones(N-1,1),1) - diag(ones(N-1,1),-1)); delta_y_M = 1/(2*Dy)*(diag(ones(M-1,1),1) - diag(ones(M-1,1),-1)); % Cast U as a vector U_vec = U_mat(:); % Mixed derivative operator A = kron(delta_y_M,delta_x_N); U_xy_num = A*U_vec; U_xy_matrix = reshape(U_xy_num,N,M); subplot(1,2,1) contourf(X,Y,U_xy_matrix) colorbar title 'Numeric U_{xy}' subplot(1,2,2) contourf(X,Y,U_xy) colorbar title 'Analytic U_{xy}'

You can obviously create the matrix yourself, but in Matlab there is tridiag for this purpose.

For example

>> full(gallery('tridiag',5,-1,2,-1))

ans =

2 -1 0 0 0 -1 2 -1 0 0 0 -1 2 -1 0 0 0 -1 2 -1 0 0 0 -1 2

Using sparse functionality available in MATLAB to generate finite difference approximation matrix is a good option.. It saves lot (indeed very much) of memory...