The OP previously answered about that:
This was actually true at that time, as x was unused, but the OP didn’t care about cleaning and making the code more safe. And now that x is used, the bug strikes back… Once again this shows the importance of being rigorous when programing.
Also, func is called with 6 arguments, whereas EDO is defined with 7 mandatory arguments. This would have been caught by the compiler if explicit interfaces had been declared.
These are contradictory statements.
When programming it is very important to keep your work organized and clear. Also, I suspect the compiler would have pointed out many of the problems you have in your code if you would adopt the style(s)/features suggested by many others. You should learn about and use
modules, as they make it easy for the compiler to ensure you are calling procedures correctlyintent, as they add clarity to both the reader and the compiler regarding your intentions with each argumentinterface blocks for procedure dummy arguments. This one will be highly valuable in your caseIn a philosophical sense, programming is an exercise in defining in a very clear and precise way the problem and solution. There is no method of being vague to take a shortcut. You must at the very least have a very clear understanding in your own mind, but even better if the code is clear enough for future readers to gain that same understanding.
I have made a different approach for the problem.
I had a code that worked well with V=ctt, but now that V(x) I am not getting the values I am expecting. This time I have defined the potential as a subroutine. My current code is:
PROGRAM P82
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, E3, E4, a, b, coef, V, error, phi1_f, phi2_f, phi3_f, phi4_f, integral
DOUBLE PRECISION:: dx, x, L, delta, V0, phi_propia1, phi_propia2, phi_propia3, valor1, valor2, valor3, E5, E6
DOUBLE PRECISION, allocatable:: phi1(:), phi2(:),phi3(:),phi4(:)
INTEGER:: N, nequ, i
coef= 3.80995d0 !eV*A^2
L=14d0
a=-L/2d0 !A
b=L/2d0 !A
V0=-50d0 !eV
delta= 0.4d0 !A
!OBRIM EL DOCUMENT
OPEN(10, file='p82.dat')
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!APARTAT 1
write(10,*) "Apartat 1"
N=400 !Número de passos
! Creació de vectors per emmagatzemar les solucions
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, V0, delta, a, b, N, phi1_f, phi1, integral, coef)
call integrate(nequ, E2, V0, delta, a, b, N, phi2_f, phi2, integral, coef)
call integrate(nequ, E3, V0, delta, a, b, N, phi3_f, phi3, integral, coef)
call integrate(nequ, E4, V0, delta, a, b, N, phi4_f, phi4, integral, coef)
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
write(10,*) ""
write(10,*) ""
CLOSE(10)
END PROGRAM P82
!SUBRUTINE RUNGE-KUTTA D'ORDRE 4
SUBROUTINE RungeKutta4order(nequ, dx, x0, yyin, yyout, func, E, V0, delta, coef)
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 :: x0, dx, E, V, V0, delta, coef, x
call func(nequ, x0, yyin, k1, E, V0, delta, coef)!k1 = f(x0, y0)
yaux = yyin + (dx * k1 * 0.5d0)
x = x0 + 0.5d0 * dx
call func(nequ, x0, yaux, k2, E, V0, delta, coef)!k2
yaux = yyin + (dx * k2 * 0.5d0)
x = x0 + 0.5d0 * dx
call func(nequ, x0, yaux, k3, E, V0, delta, coef)!k3
yaux = yyin + (dx * k3)
x = x0 + dx
call func(nequ, x0, yaux, k4, E, V0, delta, coef) !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, V0, delta, coef)
IMPLICIT NONE
DOUBLE PRECISION:: x, yin(nequ), dyout(nequ), E, V0, delta, V, coef
INTEGER:: nequ
call potencial(x, V0, delta, V)
dyout(1)= yin(2)
dyout(2) = 2D0*(V-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, V0, delta, a, b, N, phi, phi_list, integral, coef)
IMPLICIT NONE
DOUBLE PRECISION:: E, phi, a, b, V0, delta, E0, integral, V, x, coef
DOUBLE PRECISION, DIMENSION(N):: phi_list, phi2_list
DOUBLE PRECISION:: dx, yyin(nequ), yyout(nequ)
INTEGER:: nequ, i, N
EXTERNAL:: EDO, potencial
!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 !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(1)**2
phi_list(i)=yyin(1)
call potencial(x, V0, delta, V)
call RungeKutta4order(nequ, dx, x, yyin, yyout, edo, E, V0, delta, coef)
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 potencial(x, V0, delta, V)
IMPLICIT NONE
DOUBLE PRECISION:: x, V0, delta
DOUBLE PRECISION:: V
V=V0*sinh(2d0)/(cosh(2d0)+cosh(x/delta))
END SUBROUTINE potencial
I am now getting a similar plot but I am getting values of phi much higher than expected.
This is my current figure.
The next step of the exercise is to study how many iterations I have to do on the shooting method to get to the eigenvalue E3. To do so I have programed the same subroutine tir adjusting de arguments V0 and delta.
They make me calculate this convergence with 6 energies. For some reason, using E1 and E2 the subroutine never converges. I know the shooting method isn’t the bes for this case, but why does it calculate it properly with the other energies and not with this pair?
The code I added is:
SUBROUTINE tir(nequ, E1, E2, V0, delta, error, a,b,N,E3, phi3, coef)
IMPLICIT NONE
DOUBLE PRECISION:: E1, E2, E3, a, b, V0, delta, phi_list(N), coef
DOUBLE PRECISION:: phi1, phi2, phi3, dx, x, error, yyin(nequ), yyout(nequ), integral
INTEGER:: i, nequ, N, iteR
EXTERNAL:: V
!PAS 1 i 2: Definim dues energies i integrem fins phi1(1), phi2(1)
iter=100
do i = 1, iter
call integrate(nequ, E1, V0, delta, a, b, N, phi1, phi_list, integral, coef)
!print*, E1, E2
call integrate(nequ, E2, V0, delta, a, b, N, phi2, phi_list, integral, coef)
call secant(E1, E2, phi1, phi2, E3)
call integrate(nequ, E3, V0, delta, a, b, N, phi3, phi_list, integral, coef)
!Escrivim en el document el valor de E3 per cada ietració per estudiar-ne a convergència
write(10, *) i, E3
if (abs(phi3).LT.error) then
write(*,*) 'Ha convergit', i, E3
exit
else
E1 = E2
E2 = E3
end if
end do
END SUBROUTINE tir
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
And at the main program I am calling:
write(10,*) "Apartat 2"
E5 =-4d0 !eV
E6 = -3.5d0 !eV
error = 1.0d-6
call tir(nequ, E1, E2,V0, delta, error, a,b,N,Evalor1, phi1_p, coef)
call tir(nequ, E3, E4,V0, delta, error, a,b,N,Evalor2, phi2_p, coef)
call tir(nequ, E5, E6,V0, delta, error, a,b,N,Evalor3, phi3_p, coef)
Thank you very much for all your help.
EDIT: I don’t understand why for higher energies the method doesn’t converge