How to make 2D matrix from 3D matrix to solve 2D nonlinear system?

1 visualización (últimos 30 días)
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 comentarios
Bjorn Gustavsson
Bjorn Gustavsson el 4 de Jul. de 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.
Nicola Pintus
Nicola Pintus el 4 de Jul. de 2022
Yes, the elements come from the 3D heat equation. I'm editing the question, with more specifications about the problem.

Iniciar sesión para comentar.

Respuesta aceptada

Bjorn Gustavsson
Bjorn Gustavsson el 4 de Jul. de 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

Más respuestas (1)

Torsten
Torsten el 4 de Jul. de 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:

Categorías

Más información sobre Logical en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by