Reshaping of 3-dimensional array to construct an isosurface
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
Bruno Luong
le 4 Août 2023
Modifié(e) : Bruno Luong
le 4 Août 2023
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 commentaires
Bruno Luong
le 4 Août 2023
Modifié(e) : Bruno Luong
le 4 Août 2023
Supposingly the expected figure is similar to this
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1449992/image.png)
Voir également
Catégories
En savoir plus sur Linear Algebra dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!