Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB Code: Stationary Schrodinger EQ, E Spec, Eigenvalues

  1. Dec 7, 2018 #1
    Hello everyone,

    For weeks I have been struggling with this quantum mechanics homework involving writing a code to determine the energy spectrum and eigenvalues for the stationary Schrodinger equation for the harmonic oscillator. I can't find any resources anywhere. If anyone could help me get started, get my matrices and equations set up, or has worked a similar problem/written a similar code before, any help would be greatly appreciated!


    Thanks in advance!

    1. The problem statement, all variables and given/known data
    media%2F942%2F942060f1-5b27-46f2-99b7-97012f312b07%2Fimage.png
    media%2F790%2F790771bf-5a8e-4cc4-a1d7-a160db35df91%2Fimage.png



    2. Relevant equations
    Included in image above

    3. The attempt at a solution
    Her is the code I have written so far. I'm not sure if this is even close or on the right track
    % Stationary Schrodinger Equation - QHO

    clear
    clc

    hbar = 6.58E-16;
    f = 400E-9;
    w = 2*pi*f;
    m = 1;
    N = 101;
    a = 0.1;
    n = 1:N;
    r = n*a;
    l = 1;
    x = a;
    mwhbar = m*w*hbar;
    y = r;

    % Operators (Matrices)
    T = diag(-1*ones(1,N-2),2) + diag(2*ones(1,N-1),1) + diag(-1*ones(1,N),0);
    K = (1/(2*a^2)) * T; % Kinetic Energy Matrix

    Veff = -(1./r) + l*(l + 1)./(2*(r.^2));
    V = diag(Veff);
    U = (1/(2*a^2)) * V; % Potential Energy Matrix

    % Equations
    H = -(1/2*a^2)*(eigen_f(n+1) - 2 * eigen_f(n) * eigen_f(n-1));
     
  2. jcsd
  3. Dec 7, 2018 #2
    UPDATE:

    Here is an update on the code I've been working on. It is probably closer to where I'm supposed to be headed. I just have problems when it comes to calculating V and Veff


    clear
    clc

    N = 101;
    a = .01;
    n = 1:101;
    r(1:N) = n.*a;
    l = input('Enter angular quantum number l: ');
    T = zeros(N);
    V = zeros(N);
    Veff = zeros(N);


    for i = 1:N
    for j = 1:N
    if i == j
    T(i,j) = 2/(2*a^2);
    elseif i == j-1 || j == i-1
    T(i,j) = -1/(2*a^2);
    end
    end
    end


    for i = 1:N
    Veff(i,i) = (1/r) + (l*(l+1))./(2.*r.^2);
    V(i,i) = V(i,i).*Veff;
    end

    H = T + V;

    [Psi,Energy] = eig(H);

    E = diag(Energy);

    fprintf('Lowest Eigenvalues:\n');
    disp(E(1:3))


    Psi_0 = Psi(:,1);
    Psi_0_Squared = Psi_0.^2;

    Psi_1 = Psi(:,2);
    Psi_1_Squared = Psi_1.^2;


    Range_a = a:a:N*a;


    Integral_0 = trapz(Range_a,Psi_0_Squared);
    Integral_1 = trapz(Range_a,Psi_1_Squared);


    Normalized_Psi_0 = 1/sqrt(Integral_0)*Psi_0;
    Normalized_Psi_1 = 1/sqrt(Integral_1)*Psi_1;

    subplot(3,1,1);
    plot(Range_a,Normalized_Psi_0,'r');
    title('Wave Function: Ground State Hydrogen Atom');
    xlabel('x');
    ylabel('Psi_0(x)');


    subplot(3,1,2);
    plot(Range_a,Normalized_Psi_1,'r');
    title('Wave Function: First Excited State Hydrogen Atom');
    xlabel('x');
    ylabel('Psi_1(x)');


    subplot(3,1,3);
    plot(Range_a,Normalized_Psi_1,'r');
    title('Wave Function: First Excited State Hydrogen Atom');
    xlabel('x');
    ylabel('Psi_2(x)');
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?