reshaping multi dimensional array and central difference method
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Sandeep Parameshwara
le 8 Fév 2020
Commenté : Sandeep Parameshwara
le 8 Fév 2020
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).
![pic.PNG](https://www.mathworks.com/matlabcentral/answers/uploaded_files/270520/pic.png)
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
0 commentaires
Réponse acceptée
J. Alex Lee
le 8 Fév 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
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Matrix Indexing dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!