Script has no errors, but no plot is given

4 次查看（过去 30 天）
Sergio
Sergio 2024-7-10
回答： Mathieu NOE 2024-7-10
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 个评论
Aquatris
Aquatris 2024-7-10
编辑：Aquatris 2024-7-10
You are not calling the function which you called exp() that produces the plot so why do you expect a plot?
Also never use built-in function/variable names as custom function names. exp() already exist in matlab.
Sergio
Sergio 2024-7-10
编辑：Sergio 2024-7-10
The problem lies in "A(i, i) = -2 / (dr^3) + 3 / (r(i)^2 * dr) + e^(r(i)^2) / (i * hbar^3 / (sqrt(2*m)^3));" I defined the function exp(r ) since it should be connected to the line here. Should I remove the function calling at the top, and rewrite exp(r(i)^2) instead? Did that, then I got a new error related to eigvecs which is unknown.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

KSSV
KSSV 2024-7-10
Ran in:
% 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');

更多回答（1 个）

Mathieu NOE
Mathieu NOE 2024-7-10
Ran in:
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');

类别

在 Help Center 和 File Exchange 中查找有关 General Physics 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by