CODEWITHSUNDEEP
fortranfortran90fortran95
Slow learner
Slow learner

Reputation: 11

Floating point overflow

I am getting floating point overflow error in this part of code. Can any of you guys help me to find out the reason.

do j=1,ny-1
  do i=1,nx-1

    sum = 0.0d0

    do k=0,1000
      n=2.0d0*dfloat(k)+ 1.0d0
      sum = sum + ((dsinh(n*pi*x(i))*dcos(n*pi*y(j)))/((n*n*pi*pi)*dsinh(2*n*pi)))
    end do

    ue(i,j)= (x(i)/(4.0d0))- 4.0d0*sum

  end do
end do

Upvotes: 0

Views: 3286

Answers (1)

RussF
RussF

Reputation: 721

The problem is the intermediate term dsinh(2*n*pi). Consider k=1000. Then n=2001 so we need to evaluate dsinh(2001*pi) which is about 0.5*exp(6286) or over 10^2700! This is far higher than any number that can be represented in double precision. You need to reevaluate the the way you are calculating the sum. The term dsinh(n*pi*x(i)) is problematic too.

My guess is that some sort of aysmptotic expansion is required for the robust evaluation of the quotient dsinh(n*pi*x(i))/dsinh(2*n*pi). For 0<x(i)<2 this term should behave as exp(n*pi*(x(i)-2)) as n becomes large. This is will be well behaved.

Upvotes: 5

Related Questions