forums.silverfrost.com

[ognik]

Hi - new to silverfrost and also learning FORTRAN. I have a short Fortran 95 program which hits an overflow because 2 variables get the same value at one point. Should be easy to manage that, right? But all my comparisons don't work, see code below for the latest try....

The statement that should pick it up is:
if (abs(tanh(SthisX)-tanh(SlastX)).EQ.zero) ......
but whenever tanh(SthisX) becomes the same value as tanh(SlastX), (both =1.000E0) the program crashes with an overflow - AFTER the if condition. IE the If condition doesn't pick it up, but the calculation after that does and divides by 0. I have also tried testing the difference against my 'zero' parameter, same result.

Any ideas? Thanks

--------------------------------------------------
program secTanHroots
!Summary: find positive root of tanh(C) using Secant method
implicit none
integer, parameter :: dp = selected_real_kind(15, 307)
real(kind=dp) :: guess, SthisX, SnextX, SlastX, Sdiff
real(kind=dp), parameter :: zero=TINY(0.0), infinity=HUGE(0.0), root=0.0
real :: C=5
integer :: iter
integer, parameter :: n=58
character(len=n ) :: stars=REPEAT('*',n)

!----- header ----
print *, " "
print *, 'Program secTanHroots starts'
print 5, 'Enter a guess to estimate the positive root of Tanh(', C,')'
print *, ' (enter a negative number to stop)'
read (*,*) guess
print *, stars
print *, " "
print 10, 'Iter.','X(i) ', 'X(i+1)"root"', 'tanh(SthisX)','tanh(SlastX)', 'diff'
!----- end header ----

do while (guess.GE.0) ! allows trials of diff. values of guess
!----- Initialise ----
SthisX= guess
Sdiff=1
SlastX=1
iter=0
Do while (abs(Sdiff).GT.zero)
iter=iter+1
print *, tanh(SthisX)-tanh(SlastX)
if (abs(tanh(SthisX)-tanh(SlastX)).EQ.zero) then
print *, '++++++++++++++++ overflow +++++++++++++++'
go to 100
end if
SnextX= SthisX-tanh(SthisX)*((SthisX-SlastX)/(tanh(SthisX)-tanh(SlastX))) !Secant method
Sdiff = root-SnextX
print 20, iter, SthisX, SnextX, tanh(SthisX),tanh(SlastX), Sdiff

SlastX=SthisX
SthisX=SnextX
if (iter.GT.20) then ! to prevent overrun
print *, 'forced stop, cntr > 20'
stop
end if
end do

100 print *, stars
print *, " "
print *, 'Enter another guess (or negative number to stop)'
read *, guess
print *, " "
end do

print *, 'User selected end'
print *, stars
5 format (A, F3.1,A)
10 format (A, 4X,5(A, 6X))
20 format (I2, 5(3X, ES14.6))
end program secTanHroots

[mecej4]

It appears that your program is written to solve tanh(x) = 0, which has just one root, x = 0. If you really meant to solve that trivial equation instead of the more likely equation tanh(x) = C, given a value of C not equal to zero, there was no need to write a program. Your program is working reasonably, but since the estimates of the root approach zero, the secant formula will involve calculation of 0 / 0 when x is too small to represent in a machine representation.

In any case, you do not need to write a program to obtain the root, since we know that 2 arctanh(C) = ln [ (1 + C) / (1 - C)].

[davidb]

A few ideas ...

It would be a good idea to try to indent the code to make it easier to see its structure.

Do you really want to test for the absolute difference in the two function values to be equal to your "zero" parameter? Should you replace .EQ. with .LE. ?

Your constants are default real kind and not real(kind=dp). You need to use the following to get the precision you want:

[ognik]

Thanks for replies, all very useful.
. My code is indented in silverfrost, but it gets lost in the copy&paste to the forum - how do I keep that formatting please?
. It is a trivial problem - designed only to demonstrate the method using a simple function, for the 'single root'. I'm learning Fortran as part of a degree by correspondence, so the extra support is really helpful.
. Testing with .LT. instead of .EQ. did the trick , I will take comparing reals much more seriously now
. Thanks for the new knowledge of '_dp', the course I am doing kind of assumes we can teach ourselves fortran so all pure fortran advice is much appreciated.

[ognik]

Hi again, as an aside - I was thinking about dp, does anyone know why we declare it like this? integer, parameter :: dp = selected_real_kind(15, 307)
Why use 'integer'? I think because I don't understand it well enough, it looks a little clumsy?

[mecej4]

Why is dp declared integer? Because kind numbers are integers that serve to define the size of the type. The particular selection of numbers depends on the compiler. Some compilers use 4 and 8 for the kind number of 4-byte and 8-byte reals. FTN95 uses kind=1 and 2 for these types. The name 'dp' is user selected, as it is suggestive of 'double precision'.

[IanLambley]

To retain the indents, use the "code" button before and after the pasted info. You will soon find that there is a limit to the length of a code section and you may have to post it as several individual posts.
Ian

[davidb]

[ognik]

A belated thanks to all, I'm making progress