Fortran Help with fortran numerical precision

  • Thread starter Cayo
  • Start date
3
0
Hello. I need help to understand why my code is not giving me double numerical precision.
In this piece of code, I create a function that is analytically normalized, but when I calculate the numerical normalization factor, it seems to be with single precision. Here it is:

Fortran:
program dyna
implicit none
integer i
integer, parameter :: Nq1=61
real(kind=8), parameter :: saw2au=1822.889950851334d0 
real(kind=8), parameter :: nitro=14.0067d0*saw2au 
real(kind=8), parameter :: mredu=nitro**2.0d0/(2.0d0*nitro)
real(kind=8) :: e0,pi,ch,x,x0,stepX,w,expo,c0,c1
complex(kind=8) :: soma,vec0(Nq1)
pi = 3.141592653589793d0
e0=0.005d0
x0=2.09970623d0 
stepX=0.01d0
w=2.d0*e0
c0 = ((mredu*w)/pi)**(1.d0/4.d0)
c1 = (4.d0/pi*(mredu*w)**3.d0)**(1.d0/4.d0)
do i=1,Nq1
  ch=(i-(Nq1-1.d0)/2.d0-1.d0)*stepX
  x=x0+ch 
  expo = dexp(-(x-x0)**2.d0*mredu*w/2.d0)
  vec0(i) = c0 * expo
end do
!----------------------------------------!
!normalizing                             !
soma=0.0d0                               !
do i=1,Nq1                               !
  soma=soma+conjg(vec0(i))*vec0(i)*stepX !
end do                                   !
vec0=vec0/sqrt(soma)                     !
!----------------------------------------!
write(*,'(a24,2(f23.15))')'normalization constant =',soma
end program
I get 0.999998932390816, and it should be 1.0d0.
I use complex vectors because later in the code they will become complex.
I don't understand where is the problem. Can someone help me, please?
Thanks in advance,
Cayo
 

DrClaude

Mentor
6,867
3,007
kind=8 doesn't guarantee that you will get double precision. The Fortran 90 standard doesn't specify how kind is to be implemented, so it is compiler-dependent.

If you want to insure that all variables are declared as double precision, you can use
Fortran:
integer, parameter :: dp = kind(1.d0)

real(kind=dp) :: x
complex(kind=dp) :: z
 
3
0
Thank you for the reply, you are correct.
But, I changed what you suggested and it didn't changed the result. Any other idea?
 

DrClaude

Mentor
6,867
3,007
I looked at your code in more details, and I don't understand why you say that
it should be 1.0d0.
You are approximating an integral from ##-\infty## to ##\infty## using a finite grid. Why would you expect that to be an exact? Check the values of vec0(1) and vec0(Nq1).
 
3
0
Hello DrClaude,
You are right. Thats the problem, the size of the grid is too small. Thank you for the help.
 

Want to reply to this thread?

"Help with fortran numerical precision" You must log in or register to reply here.

Physics Forums Values

We Value Quality
• Topics based on mainstream science
• Proper English grammar and spelling
We Value Civility
• Positive and compassionate attitudes
• Patience while debating
We Value Productivity
• Disciplined to remain on-topic
• Recognition of own weaknesses
• Solo and co-op problem solving
Top