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 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!