How to make 2D matrix from 3D matrix to solve 2D nonlinear system?
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
Bjorn Gustavsson
2022-7-4
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.
采纳的回答
Bjorn Gustavsson
2022-7-4
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
更多回答(1 个)
Torsten
2022-7-4
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:
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!