Script has no errors, but no plot is given
23 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Sergio
le 10 Juil 2024 à 15:45
Réponse apportée : Mathieu NOE
le 10 Juil 2024 à 16:34
Hi , I wonder what is wrong here with the script. Nothinig happens when it should plot something... Can someone help? Thanks
% Constants
hbar = 1.0545718e-34; % Planck's constant
m = 9.10938356e-31; % Electron mass
l = 1; % angular momentum quantum number
E_ion = 1; % Ionization energy, example value
%Functions
function e =exp(r)
% Discretization parameters
r_min = 0.01; % Avoid division by zero
r_max = 10;
N = 100; % Number of grid points
r = linspace(r_min, r_max, N);
dr = r(2) - r(1);
% Initialize the R(r) vector
R = zeros(N, 1);
% Set up the finite difference matrix
A = zeros(N, N);
for i = 3:N-2
A(i, i-2) = -1 / (2 * dr^3);
A(i, i-1) = 2 / (dr^3) - 3 / (r(i) * dr^2);
A(i, i) = -2 / (dr^3) + 3 / (r(i)^2 * dr) + e^(r(i)^2) / (i * hbar^3 / (sqrt(2*m)^3));
A(i, i+1) = -2 / (dr^3) + 3 / (r(i) * dr^2);
A(i, i+2) = 1 / (2 * dr^3);
end
% Apply boundary conditions (example: R(0) = 0 and R(r_max) = 0)
A(1,1) = 1;
A(N,N) = 1;
% Solve the eigenvalue problem
[~, D] = eig(A);
% The eigenvalues are the diagonal elements of D
eigenvalues = diag(D);
% The solution R(r) corresponds to the eigenvector with eigenvalue closest to E_ion
[~, idx] = min(abs(eigenvalues - E_ion));
R = eigvecs(:, idx);
% Plot the solution
plot(r, R);
xlabel('r');
ylabel('R(r)');
title('Radial Wavefunction');
end
2 commentaires
Réponse acceptée
KSSV
le 10 Juil 2024 à 16:26
% Constants
hbar = 1.0545718e-34; % Planck's constant
m = 9.10938356e-31; % Electron mass
l = 1; % angular momentum quantum number
E_ion = 1; % Ionization energy, example value
%Functions
% Discretization parameters
r_min = 0.01; % Avoid division by zero
r_max = 10;
N = 100; % Number of grid points
r = linspace(r_min, r_max, N);
dr = r(2) - r(1);
% Initialize the R(r) vector
R = zeros(N, 1);
% Set up the finite difference matrix
A = zeros(N, N);
for i = 3:N-2
A(i, i-2) = -1 / (2 * dr^3);
A(i, i-1) = 2 / (dr^3) - 3 / (r(i) * dr^2);
A(i, i) = -2 / (dr^3) + 3 / (r(i)^2 * dr) + exp(r(i)^2) / (i * hbar^3 / (sqrt(2*m)^3));
A(i, i+1) = -2 / (dr^3) + 3 / (r(i) * dr^2);
A(i, i+2) = 1 / (2 * dr^3);
end
% Apply boundary conditions (example: R(0) = 0 and R(r_max) = 0)
A(1,1) = 1;
A(N,N) = 1;
% Solve the eigenvalue problem
[eigvecs, D] = eig(A);
% The eigenvalues are the diagonal elements of D
eigenvalues = diag(D);
% The solution R(r) corresponds to the eigenvector with eigenvalue closest to E_ion
[~, idx] = min(abs(eigenvalues - E_ion));
R = eigvecs(:, idx);
% Plot the solution
plot(r, R);
xlabel('r');
ylabel('R(r)');
title('Radial Wavefunction');
0 commentaires
Plus de réponses (1)
Mathieu NOE
le 10 Juil 2024 à 16:34
again...
let's fix it
I cannot check here if e^(r(i)^2) (probably wrong) must be replaced by exp(r(i)^2) - I let you double check that point , but that is actually what I did
then, you have to change one more line to define eigvecs
[~, D] = eig(A); => [eigvecs, D] = eig(A);
% Constants
hbar = 1.0545718e-34; % Planck's constant
m = 9.10938356e-31; % Electron mass
l = 1; % angular momentum quantum number
E_ion = 1; % Ionization energy, example value
% Discretization parameters
r_min = 0.01; % Avoid division by zero
r_max = 10;
N = 100; % Number of grid points
r = linspace(r_min, r_max, N);
dr = r(2) - r(1);
% Initialize the R(r) vector
R = zeros(N, 1);
% Set up the finite difference matrix
A = zeros(N, N);
for i = 3:N-2
A(i, i-2) = -1 / (2 * dr^3);
A(i, i-1) = 2 / (dr^3) - 3 / (r(i) * dr^2);
A(i, i) = -2 / (dr^3) + 3 / (r(i)^2 * dr) + exp(r(i)^2) / (i * hbar^3 / (sqrt(2*m)^3));
A(i, i+1) = -2 / (dr^3) + 3 / (r(i) * dr^2);
A(i, i+2) = 1 / (2 * dr^3);
end
% Apply boundary conditions (example: R(0) = 0 and R(r_max) = 0)
A(1,1) = 1;
A(N,N) = 1;
% Solve the eigenvalue problem
[eigvecs, D] = eig(A);
% The eigenvalues are the diagonal elements of D
eigenvalues = diag(D);
% The solution R(r) corresponds to the eigenvector with eigenvalue closest to E_ion
[~, idx] = min(abs(eigenvalues - E_ion));
R = eigvecs(:, idx);
% Plot the solution
plot(r, R);
xlabel('r');
ylabel('R(r)');
title('Radial Wavefunction');
0 commentaires
Voir également
Catégories
En savoir plus sur Annotations 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!