Reshaping of 3-dimensional array to construct an isosurface
조회 수: 1 (최근 30일)
이전 댓글 표시
I'm having an issue finalizing the code that I've ran past line 40 (code below). After reshaping my matrices to perform operations, I am unable to reshape back into an appropriate size to create an isosurface. I'm thinking I may need to slice, or perhaps I'm glossing over something entirely simple. I'm also having issues running the eig function under GPU execution for computational efficiency, but that's tertiary to the problem at hand. If you wish to execute the code to get an idea, I would suggest lowering the N value, as I have, to below 50.
Code:
%establishing iteration size
N = 25;
%establishing meshgrid spacing of appropriate bohr radius to avoid infinite
%square well
x = linspace(-25,25,N);
y = linspace(-25,25,N);
z = linspace(-25,25,N);
[X,Y,Z] = meshgrid(x,y,z);
dx = diff(x(1:2));
%estimating potential of hydrogen atom
e = 1*10^(-10);
d = (sqrt(X.^2 + Y.^2 + Z.^2 + e)).^-1;
n = -dx^2;
V = d.*n;
%second partial using kronecker product
D = ones(N,1);
D = [D, -2*D, D];
A = spdiags(D,-1:1,N,N);
%A = full(A);
SPA = sparse(kron(A,A));
%kinetic operator
T = -.5*sparse(kron(SPA,A));
V = reshape(V,1,N^3);
U = sparse(diag(V));
%hamiltonian
H = T + U;
%%H = gpuArray(H);
%calculating discrete eigenstates(energies) in three dimensions, assuming
%no progression in time
[L,C,R] = eigs(H,25,'smallestreal');
%where I'm running into issues
L = reshape(L,1,N^3);
isosurface(x,y,z,L);
댓글 수: 0
답변 (1개)
Bruno Luong
2023년 8월 4일
편집: Bruno Luong
2023년 8월 4일
I just fix the code without knowing what is really your aim. I'm surprise you do all kind of kron and eigs and get lost in the underlined dimensions of the result.
I also remove the unnecessary conversion back and forth between full ans sparse.
%establishing iteration size
N = 25;
%establishing meshgrid spacing of appropriate bohr radius to avoid infinite
%square well
x = linspace(-25,25,N);
y = linspace(-25,25,N);
z = linspace(-25,25,N);
[X,Y,Z] = meshgrid(x,y,z);
dx = diff(x(1:2));
%estimating potential of hydrogen atom
e = 1*10^(-10);
d = (sqrt(X.^2 + Y.^2 + Z.^2 + e)).^-1;
n = -dx^2;
V = d.*n;
%second partial using kronecker product
D = ones(N,1);
D = [D, -2*D, D];
A = spdiags(D,-1:1,N,N);
%A = full(A);
SPA = kron(A,A);
%kinetic operator
T = -.5*kron(SPA,A);
V = reshape(V,1,N^3);
U = diag(V);
%hamiltonian
H = T + U;
%%H = gpuArray(H);
%calculating discrete eigenstates(energies) in three dimensions, assuming
%no progression in time
NbEig = 4; % 25 for you
[L,C,R] = eigs(H,NbEig,'smallestreal');
for k=1:NbEig
subplot(2,2,k)
Lk = reshape(L(:,k),N,N,N);
isosurface(x,y,z,Lk);
end
댓글 수: 2
Bruno Luong
2023년 8월 4일
편집: Bruno Luong
2023년 8월 4일
Supposingly the expected figure is similar to this
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1449992/image.png)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!