I gets difficult to help you without understanding more the equations you are solving… Particularly what you are doing in tir() is definitely unclear to me.
I am asked to find the eigenvalues and eigenfunctions of a particle using the Schrodinger equation. To do so I created a RungeKutta4order subroutine that integrates teh schrodinger equation edo. To find the excited states I made a shooting methd tir: using two near values of E (e1 and e2) I integrate the schrodinger equation with initial conditions. I then find a new E3 which is going to give me the energy on the excited state.
Finally I am asked to fins the 4 first eigenvalues and eigenfunctions.
My current issue is found in teh fact that when I want to do a loo such as :
do i = 1, 4
E1= (pi**2/2d0 +V)
E2= ((2d0*pi)**2/2d0 +V)
print*, E1, E2
call tir(nequ, E1, E2, error, a,b,N, yyin)
enddo
I have nans or simply, it doesn’t give me the value of the convergence of tir i asked to print out. Why is that so?
How would I find the eigenvalues and eigenvectors?
I’m going to post my code straight out of what I have right now:
PROGRAM PREPA8
IMPLICIT NONE
DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b, x, T, w_n, phi
INTEGER:: N, i, iter, k
INTEGER, PARAMETER:: nequ=2
DOUBLE PRECISION:: yyin(nequ), yyout(nequ), eigenvalues(4)
EXTERNAL:: EDO, pendolSimple
OPEN(10, file='prepa8.dat')
!Dades ctts::
V=-2.4d0
pi=4d0*atan(1d0)
yyin(1)=0.02d0
yyin(2)=0.0d0
error=1.0d-5
N=500
a=0.0d0
b=1.0d0
!call EDO(nequ, x, yyin, yyout,E, V)
!print*, yyout
!do i = 1, 10
!call integrate(nequ, E, a, b, N, yyin, phi)
!print*, yyout
!yyin=yyout
!enddo
do k = 1, 4
E1= (k*pi**2/2d0 +V)
E2= (((k+1)*pi)**2/2d0 +V)
print*, E1, E2
call tir(nequ, E1, E2, error, a,b,N, yyin, E3)
eigenvalues(k)=E3
enddo
END PROGRAM PRepa8
!!!PROVA USANT PÈNDOL SIMPLE
!N=1300
!w_n = sqrt(grav/longitud)
!T= 2d0*pi/w_n
!grav = 3.71d0 !m/s^2
!longitud = 0.45d0 !m
!dx=6*T/n
!yy0(1)=0.02d0
!yy0(2)=0.0d0
!DO i = 1, N
! CALL RungeKutta4order(nequ, dx, yy0, yyout, pendolSimple)
! write(*,*) dx*i, yyout
! yy0=yyout
!ENDDO
CLOSE(10)
END PROGRAM PREPA8
!SUBRUTINE RUNGE-KUTTA D'ORDRE 4
SUBROUTINE RungeKutta4order(nequ, dx, yyin, yyout, func, E, V)
IMPLICIT NONE
INTEGER, INTENT(IN) :: nequ
DOUBLE PRECISION, DIMENSION(nequ), INTENT(IN) :: yyin
DOUBLE PRECISION, DIMENSION(nequ), INTENT(OUT) :: yyout
DOUBLE PRECISION, DIMENSION(nequ) :: k1, k2, k3, k4, yaux
DOUBLE PRECISION :: x, dx, E, V
call func(nequ, x, yyin, k1, E, V)!k1 = f(x0, y0)
!print*, k1
yaux = yyin + (dx * k1 / 2.0d0)
x = x + 0.5d0 * dx
call func(nequ, x, yaux, k2, E, V)!k2
yaux = yyin + (dx * k2 / 2.0d0)
x = x + 0.5d0 * dx
call func(nequ, x, yaux, k3, E, V)!k3
yaux = yyin + (dx * k3)
x = x + 0.5d0 * dx
call func(nequ, x + dx, yaux, k4, E, V) !k4
yyout = yyin + (dx / 6.0d0) * (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4)
!print*, yyout
END SUBROUTINE RungeKutta4order
!L'edos en aquest cas és equivalent a la del pèndol simple en aproximacions petites.
!yin=(phi, dphi/dx)
!dyout=(dphi/dx, d^2phi/dx^2)
SUBROUTINE EDO(nequ, x, yin, dyout, E, V)
IMPLICIT NONE
DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V
INTEGER:: nequ
dyout(1)= yin(2)
dyout(2) = 2D0*(V-E)* yin(1) !Correspon a la segona derivada que ens dona l'equació de Schrodinger
END SUBROUTINE EDO
SUBROUTINE integrate(nequ, E, V, a, b, N, yyin, phi)
IMPLICIT NONE
DOUBLE PRECISION:: E, phi, a, b, V
DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ)
INTEGER:: nequ, i, N
EXTERNAL:: EDO
dx=(b-a)/dble(N)
DO i = 1, N
call RungeKutta4order(nequ, dx, yyin, yyout, edo, E, V)
!write(*,*) dx*i, yyout
yyin=yyout
ENDDO
phi=yyin(1)
END SUBROUTINE integrate
SUBROUTINE secant(E1, E2, phi1, phi2, E3)
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, phi1, phi2, E3, E, V
E3=(E1 * phi2 - E2 * phi1) / (phi2 - phi1)
!print*, E3
END SUBROUTINE secant
SUBROUTINE tir(nequ, E1, E2, error, a,b,N, yyin)
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, E3, a, b, V
DOUBLE PRECISION:: phi1, phi2, phi3, dx, x, error, yyin(nequ), yyout(nequ)
INTEGER:: i, nequ, N, iter
i=0 !incialitzo
iter=500
!Pas 1 i 2
print*, E1, E2
do i = 1, iter
call integrate(nequ, E1, V, a, b, N,yyin, phi1)
!print*, E1, E2
!print*, yyin, phi1
call integrate(nequ, E2,V, a, b, N, yyin, phi2)
!print*, phi1, phi2
call secant(E1, E2, phi1, phi2, E3)
!print*, E3
call integrate(nequ, E3, V, a, b, N,yyin, phi3)
if (abs(phi3).LT.error) then
write(*,*) 'Ha convergit', i, E3
exit
else
E1 = E2
E2 = E3
end if
end do
END SUBROUTINE tir
So for the first two energies (k=1) i have good results, but for k>1 i don’t. the subroutine itself works. How can I fix it??
Please help!!!
Is the goal of the exercise to find the eigenvalues of Schrödinger equation and to learn Fortran in the process?
Both MATLAB (bvp4c) and Scipy (solve_bvp) have solvers for boundary-value problems. Perhaps you could use them as a second validation of your Fortran solution. By the looks at their examples, they should be easy to get running.
A few classic BVP solvers in Fortran can also be found on Netlib (COLSYS, COLMOD, COLNEW, MIRKDC, ACDCC). I’m not sure if there are any newer ones than BVP_SOLVER (from 2012).
But okay, your code is meant to solve for the eigenfunction of
\frac{\partial^2 \psi(r)}{\partial r^2} + \frac{2m}{\hbar^2}[E - V]\psi(r) = 0, \quad r \in [0,1]
subject to the boundary conditions \psi(0) = 0.02, \psi'(0) = 0, \psi(1) = 0, if I have inferred this correctly?
To do so, you introduced the variables \boldsymbol{y} = (y_1, y_2) = (\psi, \psi'), and rewrote the second order PDE as a system of linear first-order ODEs
\begin{align}
\frac{dy_1}{dr} &= y_2 \\
\frac{dy_2}{dr} &= 2[V - E]y_1
\end{align}
where V and E are understood to be non-dimensionalized.. You solve the forward problem with RK4, and use the secant method to iterate on the value E which meets the boundary condition at r = 1. Did I get that correctly?
2 Likes
I was talking about the decimal separator dot: all three of these are valid Fortran, all mean the same thing, and I just used the shortest: 1d-5 1.d-5 1.0d-5
1 Like
Using the equations below in MATLAB I get the following graph:
%-------------------------------------------
function dydx = mat4ode(x,y,E,V) % equation being solved
dydx = [y(2)
2*(V - E)*y(1)];
end
%-------------------------------------------
function res = mat4bc(ya,yb,E) % boundary conditions
res = [ya(2)
yb(1)
ya(1)-0.02];
end
%-------------------------------------------
function yinit = mat4init(x,k) % initial guess function
a = pi*(k - 1./2.);
yinit = [0.02*cos(a*x)
-0.02*a*sin(a*x)];
end
2 Likes
The problem (at least one problem) is in how you are updating E1 and E2, which is wrong. There are multiple cases to take into account (depending on the signs of phi1/phi2/phi3), and the update has to be different for each case.
Since this is an algorithmic issue and not a language issue, and since this is a homework, I let you think about the right solution 
… It’s not very complicated.
1 Like
But if I do separately k=1 and compute, then erase the energies with k=1 and put k=2 and compute I have the results for k=2 correctly, so the subroutine does work for k=2 and so on, the issue is when I want to update the energies on a loop and i suppose it has to do on how the energies are being passed through the subroutine inside the shooting method tir. I have thought about it and found no solution, I suppose it’s something very simple but I have no notion of what could be changed so even if it’s not complicated, I would appreciate to know what could be done…
I’ve told you that the way you are updating E1 and E2 in the tir() routine looks wrong to me, have you tried looking at that ?
Moreover the last code you have posted won’t compile (there are 2 END PROGRAM) and contains a call error (you are calling tir() with an E3 argument, which is not in the argument lists of tir())
And by continously changing and randonly changing your code you are re-introducing mistakes that you had corrected before. Notably, the integrate routine should not modify the yyin argument, because if it does then yyin is wrong for all the next calls to integrate. You have been told previsouly to add intent() to all your routine arguments to catch such mistakes, why didn’t you do so ?
I have and I don’t see how I could update E1 and E2 any other way. I do see that yyin is being modified over and over again. I tried putting it in the integrate equation instead of a parameter for the subroutine tir but I still have fault values.
Current code:
PROGRAM PREPA8
IMPLICIT NONE
DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b, x, T, w_n, phi
INTEGER:: N, i, iter, k
INTEGER, PARAMETER:: nequ=2
DOUBLE PRECISION:: yyin(nequ), yyout(nequ), eigenvalues(4)
EXTERNAL:: EDO, pendolSimple
OPEN(10, file='prepa8.dat')
!Dades ctts::
V=-2.4d0
pi=4d0*atan(1d0)
error=1.0d-5
N=500
a=0.0d0
b=1.0d0
!call EDO(nequ, x, yyin, yyout,E, V)
!print*, yyout
!do i = 1, 10
!call integrate(nequ, E, a, b, N, yyin, phi)
!print*, yyout
!yyin=yyout
!enddo
do k = 1, 4
E1= (k*pi**2/2d0 +V)
E2= (((k+1)*pi)**2/2d0 +V)
!print*, E1, E2
call tir(nequ, E1, E2, error, a,b,N,E3)
eigenvalues(k)=E3
enddo
END PROGRAM PREPA8
SUBROUTINE integrate(nequ, E0, V, a, b, N, phi)
IMPLICIT NONE
DOUBLE PRECISION:: E, phi, a, b, V, E0
DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ)
INTEGER:: nequ, i, N
EXTERNAL:: EDO
yyin(1)=0.02d0
yyin(2)=0.0d0
E=E0
dx=(b-a)/dble(N)
DO i = 1, N
call RungeKutta4order(nequ, dx, yyin, yyout, edo, E, V)
!write(*,*) dx*i, yyout
yyin=yyout
ENDDO
phi=yyin(1)
END SUBROUTINE integrate
SUBROUTINE secant(E1, E2, phi1, phi2, E3)
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, phi1, phi2, E3, E, V
E3=(E1 * phi2 - E2 * phi1) / (phi2 - phi1)
!print*, E3
END SUBROUTINE secant
SUBROUTINE tir(nequ, E1, E2, error, a,b,N,E3)
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, E3, a, b, V
DOUBLE PRECISION:: phi1, phi2, phi3, dx, x, error, yyin(nequ), yyout(nequ)
INTEGER:: i, nequ, N, iter
i=0 !incialitzo
iter=500
!print*, yyin
!print*, E1, E2
do i = 1, iter
call integrate(nequ, E1, V, a, b, N,phi1)
!print*, E1, E2
!print*, yyin, phi1
call integrate(nequ, E2,V, a, b, N,phi2)
!print*, phi1, phi2
call secant(E1, E2, phi1, phi2, E3)
!print*, E3
call integrate(nequ, E3, V, a, b, N, phi3)
if (abs(phi3).LT.error) then
write(*,*) 'Ha convergit', i, E3
exit
else
E1 = E2
E2 = E3
end if
end do
END SUBROUTINE tir
!SUBRUTINE RUNGE-KUTTA D'ORDRE 4
SUBROUTINE RungeKutta4order(nequ, dx, yyin, yyout, func, E, V)
IMPLICIT NONE
INTEGER, INTENT(IN) :: nequ
DOUBLE PRECISION, DIMENSION(nequ), INTENT(IN) :: yyin
DOUBLE PRECISION, DIMENSION(nequ), INTENT(OUT) :: yyout
DOUBLE PRECISION, DIMENSION(nequ) :: k1, k2, k3, k4, yaux
DOUBLE PRECISION :: x, dx, E, V
call func(nequ, x, yyin, k1, E, V)!k1 = f(x0, y0)
!print*, k1
yaux = yyin + (dx * k1 / 2.0d0)
x = x + 0.5d0 * dx
call func(nequ, x, yaux, k2, E, V)!k2
yaux = yyin + (dx * k2 / 2.0d0)
x = x + 0.5d0 * dx
call func(nequ, x, yaux, k3, E, V)!k3
yaux = yyin + (dx * k3)
x = x + 0.5d0 * dx
call func(nequ, x + dx, yaux, k4, E, V) !k4
yyout = yyin + (dx / 6.0d0) * (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4)
!print*, yyout
END SUBROUTINE RungeKutta4order
!L'edos en aquest cas és equivalent a la del pèndol simple en aproximacions petites.
!yin=(phi, dphi/dx)
!dyout=(dphi/dx, d^2phi/dx^2)
SUBROUTINE EDO(nequ, x, yin, dyout, E, V)
IMPLICIT NONE
DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V
INTEGER:: nequ
dyout(1)= yin(2)
dyout(2) = 2D0*(V-E)* yin(1) !Correspon a la segona derivada que ens dona l'equació de Schrodinger
END SUBROUTINE EDO
- Any declaration of NEQU as an INTEGER should precede its use in a different declaration statement as an array bound/extent.
- Variable V in TIR is not being assigned a value but is passed to INTEGRATE, and then to RUNGEKUTTA4ORDER, and then to EDO, still without an assigned value, where it is referenced, resulting in garbage.
My mistake… I was thinking about another method than the secant. That said, you know that the convergence is not guaranteed with the secant method?
For the rest, you have been reported with many errors to correct in your code and many advises about features to use to be better able to track bugs, and you don’t take them into account. So I can’t see how we can progress…
All your statements where correct, I am now having trouble normalizing and plotting the correct values of the eigenfunction, probably because the energies aren’t properly calculated. Any ideas?
Plot is supposed to look like this
In your first gnuplot graph, you are not plotting the negative part of the y-axis.
Your second gnuplot graphs look different from what your initial conditions in the integrate function are prescribing. It’s hard to judge whether \psi is zero or 0.02, but the slope, \psi' or yyin(2), doesn’t look close to zero at all. In other words, the graph you are showing us, and your Fortran program, seem to be inconsistent.
One has to be precise when coding, in the first to specify exactly what you want the computer to do, and second to assure correctness of the program.
As @PierU has mentioned, the secant method carries no guarantees of convergence.
! Secant method loop
do i = 1, iter
! evaluate phi(E1) and phi(E2)
call integrate(nequ,E1,V,a,b,N,phi_list,phi1,integral)
call integrate(nequ,E2,V,a,b,N,phi_list,phi2,integral)
! E3 = E2 - phi(E2) (E2 - E1)/(phi(E2) - phi(E1))
call secant(E1, E2, phi1, phi2, E3)
! evaluate phi(E3)
call integrate(nequ, E3, V, a, b, N,phi_list, phi3, integral)
! convergence check (absolute function value)
if (abs(phi3).LT.error) then
write(*,*) 'Ha convergit', i, E3
exit
end if
! swap values for next iteration
E1 = E2
E2 = E3
end do
Writing “imperfect” code to quickly get a result is OK, we all do that. But when facing some issues such as incorrect/unexpected results, improving the code can also help tracking and identifying the problem. With the intent() attribute the compiler can report that some variables are unexpectedly modified or used without being initialized, if the routines are in a module the compiler can check that the calls are correct, etc…
2 Likes
Hello, I have managed to program a code that returns the energy values (eigenvalues) of the Schrödinger equaions. The only thing I have left to do is normalize the eigenvectors of these values. To do so, I created a subroutine normalize that grabs the list of phi values of the eigenvalue E, calculates it’s integral and uses that to normalize the values.
I then want to plot the iteration with the normalized eigenvector.
My issue comes with the fact that I don’t know what’s failing, but the normalizing subroutine doesn’t work correctly. I have no negative values for phi which I don’t think is correct.
Please I would really appreciate some help.
My current code is:
PROGRAM PREPA8
IMPLICIT NONE
DOUBLE PRECISION:: V, E, dx, E1, E2, E3, error, pi, grav, longitud, a, b, x, T, w_n, phi3, integral
INTEGER:: N, i, iter, k
INTEGER, PARAMETER:: nequ=2
double precision, allocatable :: phi_list(:)
DOUBLE PRECISION:: yyin(nequ), yyout(nequ)
EXTERNAL:: EDO, pendolSimple
!OPEN(10, file='prepa8.dat')
!Dades ctts::
pi=4d0*atan(1d0)
error=1.0d-5
N=500
a=0.0d0
b=1.0d0
!E1=1d0
!E2=10d0
V=-2.4d0
!print*, yyout
!do i = 1, 10
!call integrate(nequ, E, a, b, N, yyin, phi)
!print*, yyout
!yyin=yyout
!enddo
allocate(phi_list(N))
do i = 1,4
E1 = (i*pi)**2/2.d0 + V +0.2d0
E2 = (i*pi)**2/2.d0 + V +0.3d0
call tir(nequ, E1, E2,V, error, a,b,N,E3, phi3)
call normalize(a, b, E, V, N, integral, phi_list)
enddo
!!!PROVA USANT PÈNDOL SIMPLE
!N=1300
!w_n = sqrt(grav/longitud)
!T= 2d0*pi/w_n
!grav = 3.71d0 !m/s^2
!longitud = 0.45d0 !m
!dx=6*T/n
!yy0(1)=0.02d0
!yy0(2)=0.0d0
!DO i = 1, N
! CALL RungeKutta4order(nequ, dx, yy0, yyout, pendolSimple)
! write(*,*) dx*i, yyout
! yy0=yyout
!ENDDO
!CLOSE(10)
END PROGRAM PREPA8
!SUBRUTINE RUNGE-KUTTA D'ORDRE 4
SUBROUTINE RungeKutta4order(nequ, dx, yyin, yyout, func, E, V)
IMPLICIT NONE
INTEGER, INTENT(IN) :: nequ
DOUBLE PRECISION, DIMENSION(nequ), INTENT(IN) :: yyin
DOUBLE PRECISION, DIMENSION(nequ), INTENT(OUT) :: yyout
DOUBLE PRECISION, DIMENSION(nequ) :: k1, k2, k3, k4, yaux
DOUBLE PRECISION :: x, dx, E, V
call func(nequ, x, yyin, k1, E, V)!k1 = f(x0, y0)
!print*, k1
yaux = yyin + (dx * k1 / 2.0d0)
x = x + 0.5d0 * dx
call func(nequ, x, yaux, k2, E, V)!k2
yaux = yyin + (dx * k2 / 2.0d0)
x = x + 0.5d0 * dx
call func(nequ, x, yaux, k3, E, V)!k3
yaux = yyin + (dx * k3)
x = x + 0.5d0 * dx
call func(nequ, x + dx, yaux, k4, E, V) !k4
yyout = yyin + (dx / 6.0d0) * (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4)
!print*, yyout
END SUBROUTINE RungeKutta4order
!L'edos en aquest cas és equivalent a la del pèndol simple en aproximacions petites.
!yin=(phi, dphi/dx)
!dyout=(dphi/dx, d^2phi/dx^2)
SUBROUTINE EDO(nequ, x, yin, dyout, E, V)
IMPLICIT NONE
DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V
INTEGER:: nequ
dyout(1)= yin(2)
dyout(2) = 2D0*(V-E)* yin(1) !Correspon a la segona derivada que ens dona l'equació de Schrodinger
END SUBROUTINE EDO
!SUBRUTINA QUE INTEGRA L'EQUACIÓ DE SHCRÖDINGER PER UNA ENERGIA DONADA
SUBROUTINE integrate(nequ, E0, V, a, b, N, phi, phi_list, integral)
IMPLICIT NONE
DOUBLE PRECISION:: E, phi, a, b, V, E0, integral
DOUBLE PRECISION, DIMENSION(N):: phi_list, phi2_list
DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ)
INTEGER:: nequ, i, N
EXTERNAL:: EDO
!PAS 1 i 2: Defineixo dues energies E1 i E2. Integrem l'equació usant el mètode de Runge-Kutta4
!Integrem tenint en compte les condicions inicials:
yyin(1)=0.02d0
yyin(2)=0.0d0
E=E0
dx=(b-a)/dble(N)
DO i = 1, N
phi2_list(i)=yyin(2)**2
phi_list(i)=yyin(2)
call RungeKutta4order(nequ, dx, yyin, yyout, edo, E, V)
!write(*,*) dx*i, yyout
yyin=yyout
ENDDO
phi=yyin(2)
call trapezoidalrule_list(a, b, N, phi2_list, integral)
END SUBROUTINE integrate
SUBROUTINE secant(E1, E2, phi1, phi2, E3)
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, phi1, phi2, E3, E, V
E3=(E1 * phi2 - E2 * phi1) / (phi2 - phi1)
!print*, E3
END SUBROUTINE secant
SUBROUTINE tir(nequ, E1, E2,V, error, a,b,N,E3, phi3)
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, E3, a, b, V, phi_list(N)
DOUBLE PRECISION:: phi1, phi2, phi3, dx, x, error, yyin(nequ), yyout(nequ), integral
INTEGER:: i, nequ, N, iter
!PAS 1 i 2: Definim dues energies i integrem fins phi1(1), phi2(1)
iter=500
do i = 1, iter
call integrate(nequ, E1, V, a, b, N, phi1, phi_list, integral)
call integrate(nequ, E2, V, a, b, N, phi2, phi_list, integral)
call secant(E1, E2, phi1, phi2, E3)
call integrate(nequ, E3, V, a, b, N, phi3, phi_list, integral)
if (abs(phi3).LT.error) then
write(*,*) 'Ha convergit', i, E3
!write(*, *) E3 !escrivim el eigenvalue
exit
else
E1 = E2
E2 = E3
end if
end do
END SUBROUTINE tir
SUBROUTINE trapezoidalrule_list(a, b, N, func, resultat)
IMPLICIT NONE
DOUBLE PRECISION:: a,b, resultat, h, x, dx
DOUBLE PRECISION, DIMENSION(N)::func
INTEGER:: i, N
resultat = 0.d0
dx = (b-a)/N !Definim el pas amb l'interval donat
do i = 0, size(func)-1
resultat = resultat + dx*func(i+1)
end do
resultat = resultat - ((func(1)+func(N))*dx/2.d0)
END SUBROUTINE trapezoidalrule_list
SUBROUTINE normalize(a, b, E, V, N, integral, llista)
IMPLICIT NONE
DOUBLE PRECISION:: llista(N), llista_n(N), integral, a, b, E, dx, x, phi, V
INTEGER:: i, N, nequ
nequ=2
dx = (b-a)/N
call integrate(nequ, E, V,a, b, N, phi, llista, integral)
write(10, *) "E=", E, "N=", N
do i = 1, N
llista_n(i) = llista(i)/sqrt(integral)
write(10,'(3F15.12)') i*dx, llista(i), llista_n(i)
end do
write(10,*) ""
write(10,*) ""
END SUBROUTINE normalize
And when I plot I have:
Thank you!
Do you see what could be missing in the normalization subroutine?
When you plot, you should probably show also the negative part of the y-axis, in gnuplot:
plot [-2:2] "prepa8.dat" u 1:3 skip 1
2 Likes
Thank you very much.
I am trying a very similar exercise but in this case V depends on x. I have defined a function that calculates this V (with the function given) and I am trying to pass the values on the subroutines, defining this function as an external. The code works well but I am not getting the values expecting and I believe it could be because of this.
The code is:
PROGRAM P82022
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, E3, E4, a, b, coef, V, phi1_f, phi2_f, phi3_f, phi4_f, integral
DOUBLE PRECISION:: dx, x, L, delta, V0
DOUBLE PRECISION, allocatable:: phi1(:), phi2(:),phi3(:),phi4(:)
INTEGER:: N, nequ, i
EXTERNAL:: V
coef= 3.80995d0 !eV*A^2
L=14d0
!Defineixo la longitud de la caixa
a=-L/2d0 !A
b=L/2d0 !A
!OBRIM EL DOCUMENT
OPEN(10, file='p8.dat')
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!APARTAT 1
write(10,*) "Apartat 1"
N=400 !Número de passos
allocate (phi1(N), phi2(N),phi3(N),phi4(N))
!Defineixo les energies
E1=-31d0 !eV
E2=-30d0 !eV
E3=-14d0 !eV
E4=-13d0 !eV
nequ=2
call integrate(nequ, E1, V, a, b, N, phi1_f, phi1, integral)
call integrate(nequ, E2, V, a, b, N, phi2_f, phi2, integral)
call integrate(nequ, E3, V, a, b, N, phi3_f, phi3, integral)
call integrate(nequ, E4, V, a, b, N, phi4_f, phi4, integral)
dx=(b-a)/dble(N)
do i = 1, N
x= a+ dx*i
write(10, *) x, phi1(i), phi2(i), phi3(i), phi4(i)
enddo
!Tanquem el document
CLOSE(10)
END PROGRAM P82022
!SUBRUTINE RUNGE-KUTTA D'ORDRE 4
SUBROUTINE RungeKutta4order(nequ, dx, yyin, yyout, func, E, V)
IMPLICIT NONE
INTEGER, INTENT(IN) :: nequ
DOUBLE PRECISION, DIMENSION(nequ), INTENT(IN) :: yyin
DOUBLE PRECISION, DIMENSION(nequ), INTENT(OUT) :: yyout
DOUBLE PRECISION, DIMENSION(nequ) :: k1, k2, k3, k4, yaux
DOUBLE PRECISION :: x, dx, E, V
EXTERNAL:: V
call func(nequ, x, yyin, k1, E, V)!k1 = f(x0, y0)
!WRITE(*, *) ' :', x, ' yyin:', yyin, ' k1:', k1
yaux = yyin + (dx * k1 / 2.0d0)
x = x + 0.5d0 * dx
call func(nequ, x, yaux, k2, E, V)!k2
!WRITE(*, *) 'Step:', 0, ' x:', x, ' yyin:', yyin, ' k2:', k2
yaux = yyin + (dx * k2 / 2.0d0)
x = x + 0.5d0 * dx
call func(nequ, x, yaux, k3, E, V)!k3
!WRITE(*, *) 'Step:', 0, ' x:', x, ' yyin:', yyin, ' k3:', k3
yaux = yyin + (dx * k3)
x = x + 0.5d0 * dx
call func(nequ, x + dx, yaux, k4, E, V) !k4
yyout = yyin + (dx / 6.0d0) * (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4)
!WRITE(*, *) 'Step:', 0, ' x:', x, ' yyin:', yyin, ' k4:', k4
!print*, yyout
END SUBROUTINE RungeKutta4order
!L'edos en aquest cas és equivalent a la del pèndol simple en aproximacions petites.
!yin=(phi, dphi/dx)
!dyout=(dphi/dx, d^2phi/dx^2)
SUBROUTINE EDO(nequ, x, yin, dyout, E, V, coef)
IMPLICIT NONE
DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V, coef
INTEGER:: nequ
dyout(1)= yin(2)
dyout(2) = 2D0*(V(x)-E)* yin(1)/coef !Correspon a la segona derivada que ens dona l'equació de Schrodinger
END SUBROUTINE EDO
!SUBRUTINA QUE INTEGRA L'EQUACIÓ DE SHCRÖDINGER PER UNA ENERGIA DONADA
SUBROUTINE integrate(nequ, E0, V, a, b, N, phi, phi_list, integral)
IMPLICIT NONE
DOUBLE PRECISION:: E, phi, a, b, V, E0, integral, x, V0, delta
DOUBLE PRECISION, DIMENSION(N):: phi_list, phi2_list
DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ)
INTEGER:: nequ, i, N
EXTERNAL:: EDO, V
!PAS 1 i 2: Defineixo dues energies E1 i E2. Integrem l'equació usant el mètode de Runge-Kutta4
!Definim les condicions inicials
yyin(1)=0 !A^(-1/2)
yyin(2)=2.0d-6 !A^(-3/2)
E=E0
dx=(b-a)/dble(N)
DO i = 1, N
x=a+dx*i
phi2_list(i)=yyin(2)**2
phi_list(i)=yyin(2)
call RungeKutta4order(nequ, dx, yyin, yyout, edo, E, V)
!WRITE(*, *) 'Step:', i, ' x:', x, ' phi:', phi_list(i), ' V(x):', V(x)
yyin=yyout
ENDDO
phi=yyin(2)
!call trapezoidalrule_list(a, b, N, phi2_list, integral)
END SUBROUTINE integrate
DOUBLE PRECISION FUNCTION V(x)
DOUBLE PRECISION:: V0, delta, x
V0=-50d0 !eV
delta= 0.4d0 !A
V=V0*sinh(2d0)/(cosh(2d0)+cosh(x/delta))
END FUNCTION V
The plot I am getting looks like:
And I should be getting
What am I doing wrong? Please help!
Am I over-using EXTERNAL for V ?
The subroutines and main program work for a constant V, why doesn’t it work for a V(x)?
I am afraid that you have to discipline yourself a bit better. When you modify a previously working code, do not break it by making careless modifications.
In your latest subroutine RungeKutta, how is x defined (i.e., given a value before being used to evaluate an expression, passed as an input argument to another subprogram, etc.) ?
1 Like