How to make 2D matrix from 3D matrix to solve 2D nonlinear system?
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Nicola Pintus
le 4 Juil 2022
Commenté : Bjorn Gustavsson
le 5 Juil 2022
Hello,
I'm trying to create an algorithm to solve a system Ax=b, where A is a NxN matrix. The starting point is the 3D heat equation and the use of the Finite Difference Method. In the 3D domain, a cube of size 50x50x50, I have the boundary Dirichlet condition in two opposite faces, while on the other 4 faces I have Neumann condition. After common discretization, I get this 3D matrix of size 50x50x50 like this
A = zeros(50,50,50);
A(1,1:3) = [-3 4 1];
A(2,2) = -1;
A(1,1:2,2) = [-4 -1];
A(2,1:3,2) = [-1 6 -1];
A(3,2:3,2) = [-1 -1];
A(1,1,3) = 1;
A(2,2:3,3) = [-1 -1];
A(3,2:4,3) = [-1 6 -1];
A(4,3:4,3) = [-1 -1];
A(3,3:4,4) = [-1 -1];
A(4,3:5,4) = [-1 6 -1];
A(5,4:5,4) = [-1 -1];
A(4,4:5,5) = [-1 -1];
A(5,4:6,5) = [-1 6 -1];
A(6,5:6,5) = [-1 -1];
A(5,5:6,6) = [-1 -1];
A(6,5:7,6) = [-1 6 -1];
A(7,6:7,6) = [-1 -1];
and go on until the last pages which are:
A(47,47:48,48) = [-1 -1];
A(48,47:49,48) = [-1 6 -1];
A(49,48:49,48) = [-1 -1];
A(50,50,48) = 1;
A(48,48:49,49) = [-1 -1];
A(49,48:50,49) = [-1 6 -1];
A(50,49:50,49) = [-1 -4];
A(49,49,50) = -1;
A(50,48:50,50) = [1 -4 3];
N = size(A,1)*size(A,2)*size(A,3);
As you can see, the 3D matrix has a repetitive scheme in every page (especially from the page 4 to the page 46). The first 3 pages and the last 3 pages are symmetric. My issue is to create a 2D squared matrix NxN in order to solve a system using a backslah operator.
Moreover, I'd like to write matrix A in a better, faster way.
Any advice would be much appreciated. Thanks for reading.
2 commentaires
Bjorn Gustavsson
le 4 Juil 2022
Is this some kind of 3-D discrete laplace-operator? You might get better help faster if you also specify what problem you're trying to solve instead of leaving it to others to interpret what it might be.
Réponse acceptée
Bjorn Gustavsson
le 4 Juil 2022
For the interior of your grid you can rather easily calculate the Laplacian operator something like this:
NperSide = 50;
[x,y,z] = meshgrid(1:NperSide);
x(:) = 1:NperSide^3;
idxR = zeros(numel(x)*7,1);
idxC = zeros(numel(x)*7,1);
vals = zeros(numel(x)*7,1);
idxCurr = 1;
for i1 = 2:(size(x,1)-1)
for i2 = 2:(size(x,2)-1)
for i3 = 2:(size(x,3)-1)
idxR(idxCurr) = x(i1,i2,i3);
idxC(idxCurr) = x(i1,i2,i3);
vals(idxCurr) = -6;
idxCurr = idxCurr + 1;
idxR(idxCurr) = x(i1,i2,i3);
idxC(idxCurr) = x(i1-1,i2,i3);
vals(idxCurr) = 1;
idxCurr = idxCurr + 1;
idxR(idxCurr) = x(i1,i2,i3);
idxC(idxCurr) = x(i1+1,i2,i3);
vals(idxCurr) = 1;
idxCurr = idxCurr + 1;
idxR(idxCurr) = x(i1,i2,i3);
idxC(idxCurr) = x(i1,i2-1,i3);
vals(idxCurr) = 1;
idxCurr = idxCurr + 1;
idxR(idxCurr) = x(i1,i2,i3);
idxC(idxCurr) = x(i1,i2+1,i3);
vals(idxCurr) = 1;
idxCurr = idxCurr + 1;
idxR(idxCurr) = x(i1,i2,i3);
idxC(idxCurr) = x(i1,i2,i3-1);
vals(idxCurr) = 1;
idxCurr = idxCurr + 1;
idxR(idxCurr) = x(i1,i2,i3);
idxC(idxCurr) = x(i1,i2,i3+1);
vals(idxCurr) = 1;
idxCurr = idxCurr + 1;
end
end
end
MD23D = sparse(idxR(vals~=0),idxC(vals~=0),vals(vals~=0));
As Torsten notes this matrix will be rather large, and most likely "challenging" to invert. You'll also have to fix the edges such that you properly respect your boundary conditions.
HTH
2 commentaires
Plus de réponses (1)
Torsten
le 4 Juil 2022
The matrix you have to create and work with is a (sparse) matrix of size (125000 x 125000).
The 3d matrix above (whatever it may represent) is of no use for the computation.
This code might be of help:
Voir également
Catégories
En savoir plus sur Logical 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!