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

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.

[Bjorn Gustavsson]

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]

Yes, the elements come from the 3D heat equation. I'm editing the question, with more specifications about the problem.

Answer 1

[Bjorn Gustavsson]

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

[Nicola Pintus]

thanks a lot, it works for me!

[Bjorn Gustavsson]

Great, happy that it helps.

Answer 2

[Torsten]

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:

https://de.mathworks.com/matlabcentral/fileexchange/27279-laplacian-in-1d-2d-or-3d?s_tid=FX_rc3_behav

[Nicola Pintus]

the link you sent me is very useful. Thanks a lot!