reshaping multi dimensional array and central difference method

Hello members,

I have the arrays

A1=[-1.1, -1, -0.9, -0.1, 0, 0.1, 0.9, 1, 1.1];
A2=[-1.1, -1, -0.9, -0.1, 0, 0.1, 0.9, 1, 1.1];

I create a 2D grid using these arrays, like below in this image(just for presentation. I mean, I have a gridded data in workspace).

I have matrix data at all these grid points (at both dark and smaller dots). I want to apply central difference scheme for matrices at dark dots. For example, lets say my matrix is 'T'. To get T_dot at (-1,-1) I have to write

T_dot(-1,-1)= (T(-0.9,-1)-T(-1.1,-1))/0.2...  % Central difference applied in 2 direction
            + (T(-1,-0.9)-T(-1,-1.1))/0.2
        

Hence,I want derivatives of T at all these dark dots(That is for 9 points here in the image).

T_dot(0,-1) = ...
    .
    .
T_dot(1,1) = ...

I want to automate this for N dimensional grids. The array value will be same in all direction. So how can I write a flexible code for this in MATLAB? My bet would be on using ndgrid, reshape functions and using a single for loop and later reshape it back. But I am kind of stuck with writing indices in the loop.

Could you kindly help me?

Thanks a lot

Answer 1

[J. Alex Lee]

Based on your example lines, you seem to have the further simplification that your delta x and delta y are uniform (in your example equal to 0.2). It's easy to generalize (and actually more readable) where dx~=dy

The finite difference coefficients matrix is very sparse, especially as N increases, so you want to use sparse matrices and assignments. Let's assume your coefficients matrix is called A rather than T_dot, with size [N,M].

First preallocate

A = sparse(N*M,N*M);

Then, kind of in the spirit of re-shaping, you can loop through the indices of A, and get the subscripts using ind2sub(). These subscripts allow you to easily find the elements that form your "cross", that is, in the internal nodes.

Once you find the nodes of interest using the subscripts, you can convert those back into the linear indices into the "unwrapped" version of A to do your assignments.

The most efficient way to populate sparse matrices (I think) is to vectorize the assignment. So you can create a temporary array s with the number of elements equal to the number of affected nodes in the "cross", i.e., 5, and write your equations into it.

Below example uses simple diffusion as an example, with D the diffusion coefficient. So rather than your central 1st order difference, it is based on 2nd order difference; should be easy to adapt to your case.

s = nan(1,5);
for k = 1:N*M
    [i,j] = ind2sub([N,M],k);
    % order of subscripts:
    % [  i,  j]
    % [i+1,  j]
    % [i-1,  j]
    % [  i,j-1]
    % [  i,j-2]
    
    % this can only work for internal nodes so you need to have some if statements
    % see below note on boundary conditions
    
    % p will be linear the index into the grid
    % which will correspond to the column index in the Jacobian
    %
    p = sub2ind([N,M],...
    [i,i+1,i-1,i  ,i  ],...
    [j,j  ,j  ,j-1,j-2]);
    
    % p = [i,j]
    s(1) = - D*2/dx/dx  + D/dy/dy;
    
    % p = [(i+1),j]
    s(2) = + D/dx/dx;
    
    % p = [(i-1),j]
    s(3) = + D/dx/dx;
    
    % p = [i,(j-1)]
    s(4) = - D*2/dy/dy;
    
    % p = [i,(j-2)]
    s(5) = + D/dy/dy;
    
    % subscript index into sparse A
    A = A + sparse(k,p,s,N*M,N*M);
    
end

You didn't address boundary conditions in your question, but this scheme makes it pretty straightforward to deal with them as well. You can detect the corners and boundaries with if statements

for k = 1:N*M
    [i,j] = ind2sub([N,M],k);
    if i==1 && j==1 % corner
    elseif i==1 % left boundary
    elseif i==N % right boundary
    elseif j==1 % top/bottom boundary
    % ...etc
    else % internal nodes
        % then you can always do i-1 and i+1 and j-1 and j+1 without out-of-bounds errors
    end
end

[Sandeep Parameshwara]

Thanks a lot for your help