2

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:

enter image description here

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!

karakfa
  • 66,216
  • 7
  • 41
  • 56
Mike Miller
  • 253
  • 6
  • 18

3 Answers3

3

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}'

Numerical and Analytic Comparison for Mixed Finite Difference Operator

Steve
  • 1,579
  • 10
  • 23
  • Thanks, this is a nice answer. What if u(x1,x2) is on a rectangular domain? – Mike Miller Mar 15 '16 at 15:37
  • Things get a little more tricky, but it's possible. You want a different number of points in each direction yes? – Steve Mar 15 '16 at 15:46
  • Yes, I would like to generalise :) I've essentially already done the first part of your suggestion. E.g. if it's discretised into a 6x4 grid, I generate two tridiagonals, 4x4 (for u_xx) and 6x6 (for u_yy) with the same structure I used in my original post. The do u*(4x4) for u_xx and (6x6)*u for u_yy. – Mike Miller Mar 15 '16 at 15:51
  • Thanks, I'll have a go! – Mike Miller Mar 15 '16 at 15:54
  • Thanks for the edit too - great help - exactly what o was looking for :) – Mike Miller Mar 18 '16 at 19:17
  • @Mike Glad to help. Also seen some similar questions before so now have something I can link to. I'm pretty sure this is order `N+M` which is not so great – Steve Mar 20 '16 at 12:52
  • Yeah, I have found my code is quite slow. I'm solving a nonlinear monge-ampere type PDE for optimal transport (with application to images). Anything roughly over a 100x100 discretisation is running slowly. Just looking for ways to speed it up. – Mike Miller Mar 20 '16 at 12:58
  • 1
    @Mike Well the matrix will be big and relatively sparse (somewhat banded, maybe width `max(N+1,M+1)`, not symmetric), so you probably want to avoid calculating it or its inverse directly. Not sure what to suggest and certainly can't think of a quick fix. – Steve Mar 20 '16 at 13:18
0

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
karakfa
  • 66,216
  • 7
  • 41
  • 56
  • Sorry, I should have been more specific in my question. I know how to create the tridiagonal matrix in matlab, it's the structure of the matrix for the partial derivatives I wondered about. Perhaps I should have asked on comp sci? – Mike Miller Mar 15 '16 at 14:24
  • sorry, I didn't pay attention, I thought you were describing the structure of this matrix. I'll think about this. Did you meant `u_ij=...` instead of `u_xy` in the equation. – karakfa Mar 15 '16 at 14:31
  • I don't think it's possible with simple matrix multiplication since the matrix indices you depend on doesn't match either row or column indices. – karakfa Mar 15 '16 at 14:37
  • Thanks :) u_xy denotes the partials w.r.t. x & y. So at a point u_{i,j} the derivation can be approximated by the above scheme I'd written. – Mike Miller Mar 15 '16 at 14:37
  • I wondered if for instance, you turned a 3x3 matrix into a 1 x 9 row vector, and multiplied by a 9 x 9 matrix (essentially getting rid of the 2D problem). It seemed quite a clunky way of doing it though and I assumed there must be something nicer out there. – Mike Miller Mar 15 '16 at 14:40
0

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

Ashik A S
  • 1