6 views (last 30 days)

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

J. Alex Lee
on 8 Feb 2020

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

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.