Witajcie programiści.
Program, który tu przedstawiam służyć powinien do obliczania wartości własnych hamiltonianu, które zależne są od pewnego współczynnika alfa. Niestety wykres, który otrzymuję jest nieprawidłowy. Prosiłbym bardzo o sprawdzenie kodu, byłbym bardzo wdzięczny.
Dodatkowo potrzebuję małej pomocy z utworzeniem histogramu z tego wykresu, który poda mi liczbę wartości w przedziale energii 3 - 3.5, który to przedział będzie podzielony na N binów. Do histogramu potrzebna jest też krzywa uwzględniająca każdy z "binów".
Oto kod:
clear all
% STAŁE--\
m=1;
h=1;
n=6; % PARAMETR POTENCJAŁU
lambda=1.5;
N=70; % LICZBA FUNKCJI BAZOWYCH
k_0=23;
B_0=1000;
B_n=0.1;
a=0.000003; %PARAMETR POTENCJAŁU
%---------------/
%ZMIENNE ZALEŻNE-\
x_max=(2*n)/lambda ;
V=a*(x_max.^(2*n))*(exp(-lambda*x_max));
q=(B_n/B_0)^(1/(N-1)); %PARAMETR POTENCJAŁU
B=zeros(N,1);
k=zeros(N,1);
B(1)=B_0;
k(1)=k_0;
%----------------/
%BAZA B, k------\
for i=2:1:N
B(i)=B(i-1)*q;
k(i)=k_0;
end
%---------------/
norma=zeros(N,1);
S=zeros(N);
T=zeros(N);
%OBLICZANIE WARTOŚCI FUNKCJI BAZOWYCH----\
for i=1:N
for j=1:i
S(i,j)=((2*factorial(2*k_0))/(B(i)+B(j)).^((2*k_0)+1));
T(i,j)=((h.^2)/m)*(((k(i)*k(j))*(factorial(k(i)+k(j)-2))/((B(i)+B(j)).^(k(i)+k(j)-1))-((k(i)*B(j)+k(j)*B(i))*(factorial(k(i)+k(j)-1))/((B(i)*B(j)).^(k(i)+k(j)))+(B(i)*B(j)*(factorial(k(i)+k(j)))/((B(i)+B(j)).^(k(i)+k(j)+1))))));
end
norma(i)=1/sqrt(S(i,i));
%NORMALIZACJA FUNKCJI S i T----------\
for j=1:i
S(i,j)=S(i,j)*norma(i)*norma(j);
S(j,i)=S(i,j);
T(i,j)=T(i,j)*norma(i)*norma(j);
T(j,i)=T(i,j);
end %---------------------------------/
%-----------------------------------------/
end
%TWORZENIE CIĄGU GEOMETRYCZNEGO WARTOŚCI WSP. ALFA----\
zakres=100;
alfa=zeros(zakres,1);
alfa_0=0.1;
alfa_n=1;
alfa_ciag_geo=(alfa_n/alfa_0).^(1.0/(zakres-1));
alfa(1)=alfa_0;
for i=2:1:zakres
alfa(i)=alfa(i-1)*alfa_ciag_geo;
end
%------------------------------------------------------/
E=zeros(N,1);
V_alfa=zeros(N);
H=zeros(N);
S_alfa=zeros(N);
figure
hold on
%OBLICZANIE WARTOSCI POTENCJALU ZALEZNEGO OD WSP. ALFA---\
for ik=1:zakres
alf=alfa(ik);
for i=1:N
for j=1:i
V_alfa(i,j)=(2*a*(alf.^(k(i)+k(j)))*(factorial(k(i)+k(j)+2*n)))/((B(i)+B(j))*alf+lambda).^(k(i)+k(j)+(2*n)+1);
V_alfa(i,j)=V_alfa(i,j)*norma(i)*norma(j);
V_alfa(j,i)=V_alfa(i,j);
end
end
%---------------------------------------------------------/
H=alf*T+V_alfa; %OBLICZANIE HAMILTONIANU
S_alfa=S/alf;
E=eig(H,S_alfa); %ENERGIA = WARTOŚCI WŁASNE HAMILTONIANU
plot(alf,E,'.') %WYKRES ZALEŻNOŚCI ENERGII od WSP. ALFA
axis([0 alf 0 V+0.1])
end
Z góry dziękuję ze wszelką pomoc i cierpliwość.
Pozdrawiam