Info

此问题已关闭。 请重新打开它进行编辑或回答。

a code using function_handle is not woring

9 次查看（过去 30 天）
Demetris
Demetris 2024-11-13，18:15
关闭： Stephen23 2024-11-13，20:41
Ran in:
matlab
% Define the domain and grid size
L = 3; % length of the domain
N = 2; % number of grid points in each direction
x = linspace(-L/2,L/2,N);
y = linspace(-L/2,L/2,N);
z = linspace(-L/2,L/2,N);
[X,Y,Z] = meshgrid(x,y,z);
dx = x(2)-x(1);
% Define the potential function
V = -1./sqrt(X.^2+Y.^2+Z.^2); % Coulomb potential
% Define the Laplacian operator using central differences
Lap = @(f) (circshift(f,[0,0,-1])+circshift(f,[0,0,1])+circshift(f,[0,-1,0])+circshift(f,[0,1,0])+circshift(f,[-1,0,0])+circshift(f,[1,0,0])-6*f)/(dx^2)*(-0.5);
% Set up the initial wave function (1s state)
psi = exp(-sqrt(X.^2+Y.^2+Z.^2)); % Gaussian wave packet
% Normalize the wave function
psi = psi./sqrt(sum(abs(psi(:)).^2*dx^3));
% Define the Hamiltonian operator
H=Lap+diag(V(:));
Operator '+' is not supported for operands of type 'function_handle'.
% Set up the time evolution parameters
dt = 0.01; % time step
%tmax = 10; % maximum time
tmax = 0.02; % maximum time
t = 0:dt:tmax;
% Solve the time-dependent Schrodinger equation using the Crank-Nicolson method
for n = 1:length(t)
psi = (eye(N^3)-0.5i*dt*H)\(eye(N^3)+0.5i*dt*H)*psi;
psi = psi./sqrt(sum(abs(psi(:)).^2*dx^3)); % renormalize
end
% Calculate the energy levels by finding the eigenvalues of the Hamiltonian
[E,D] = eigs(H,10,'sr');
E = diag(E); % extract the eigenvalues
E = E(E<0); % take only the negative energy levels (bound states)
% Print the energy levels
disp('Energy levels:');
disp(E);

此问题已关闭。

回答（0 个）

此问题已关闭。

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by