Reshaping of 3-dimensional array to construct an isosurface
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (1)
Bruno Luong
am 4 Aug. 2023
Bearbeitet: Bruno Luong
am 4 Aug. 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 Kommentare
Bruno Luong
am 4 Aug. 2023
Bearbeitet: Bruno Luong
am 4 Aug. 2023
Supposingly the expected figure is similar to this
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1449992/image.png)
Siehe auch
Kategorien
Mehr zu Linear Algebra finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!