Problem with precision solving one math series!

[guilhermemarques]

Hey guys,

I have to program a Fortran code to solve the problem: S = 1/1 + 1/2 + ... + 1/T to get the largest S possible.
Actually my fortran program is like:

program test
implicit none
real :: T = 17000000, S = 0;
do while (T<17000000)
S = S + 1/T
T = T + 1
end do

print *,"S Value = ",S


end program test

And it just don´t solve the problem, for example, with a T value 16000000 the program solves. So i think it´s the precision of the real variable...
How can i solve this problema in Fortran?

Thanks in advance.

 

[ArkM]

You can't solve any problems with this program. See: T < 17000000 is equal to false if T = 17000000...
Besides that, change real to double precision to provide ~16 digits precision instead ~7 for real.

 

[guilhermemarques]

My mistake, in the original program was:
"do while (T<17000099)"

How can i declare one double variable in Fortran?
I tried double :: var
But it doesn´t work.

 

[xwb]

It varies depending on which compiler you're using

In F77
DOUBLE PRECISION

In Intel IVF/CVF
REAL*8

In Salford/Silverfrost
REAL*3
REAL*2 takes 8 bytes

I've also seen
REAL*10
REAL*12
REAL*16

but not all compilers take them.

 

[guilhermemarques]

Thanks xwb, im using a G77 compiler.
Actually the real*8 works.
But the maximum value for T is yet small. When i tried to reach "(T<2200000000)" for example, it doesn´t work.
Error msg:

~Fonte1.f:5:
do while (T<2200000000)
^
Integer at (^) too large

My code:

program teste
implicit none
real*8 :: T =1, S= 0
print *,"Tataluga Project"
do while (T<2100000000)
S = S + 1/T
T = T + 1
end do

print *,"S value = ",S

end program teste

Any tips? Like using another compiler, another kind of precision instead of real*8, another code struct?

Thanks.

 

[mikrom]

Sum(1/t) where t goes from 1 to infinity is known as harmonic series. This infinity series is divergent but it diverges very slowly.
You can compute the sum of a series so long, until sum_of_series - sum_of_series_previous is approximately equal zero.
However, in numerical computation instead of using absolute zero it is rather used a approximation of the zero - aka machine epsilon, which is a number which fulfills the condition: macheps + 1 = macheps (analogical to 0+1 = 0)
Machine epsilon is machine dependent.

So here is an example

[COLOR=#a020f0]module[/color] my_numeric
[COLOR=#a020f0]contains[/color]
  [COLOR=#a020f0]function[/color] machine_epsilon(start) [COLOR=#a020f0]result[/color] (macheps)
    [COLOR=#0000ff]! Machine Epsilon[/color]
[COLOR=#2e8b57][b]    real[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color]([COLOR=#2e8b57][b]in[/b][/color]) :: start
[COLOR=#2e8b57][b]    real[/b][/color] :: macheps
    macheps [COLOR=#804040][b]=[/b][/color] start
    [COLOR=#804040][b]do[/b][/color] [COLOR=#804040][b]while[/b][/color] (macheps [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1.0[/color] [COLOR=#804040][b]>[/b][/color] [COLOR=#ff00ff]1.0[/color])
      macheps [COLOR=#804040][b]=[/b][/color] macheps [COLOR=#804040][b]/[/b][/color] [COLOR=#ff00ff]2.0[/color]
      [COLOR=#0000ff]!write (*,*) 'macheps =', macheps[/color]
    [COLOR=#804040][b]end do[/b][/color]   
  [COLOR=#a020f0]end function[/color] machine_epsilon 

  [COLOR=#a020f0]subroutine[/color] harmonic_series(macheps)
    [COLOR=#0000ff]! t and sum(1/t) aka harmonic series[/color]
[COLOR=#2e8b57][b]    real[/b][/color] :: macheps, s, s_old, t    
    s_old [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0.0[/color]
    s [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1.0[/color]
    t[COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#0000ff]!write (*,*) 't = ', t,':         s = ', s, ' s_old = ', s_old [/color]
    [COLOR=#804040][b]do[/b][/color] [COLOR=#804040][b]while[/b][/color] ((s [COLOR=#804040][b]-[/b][/color] s_old) [COLOR=#804040][b]>[/b][/color] macheps)
      t [COLOR=#804040][b]=[/b][/color] t [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
      s_old [COLOR=#804040][b]=[/b][/color] s
      s [COLOR=#804040][b]=[/b][/color] s [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1.0[/color][COLOR=#804040][b]/[/b][/color]t
      [COLOR=#0000ff]!write (*,*) 't = ', t,':         s = ', s, ' s_old = ', s_old [/color]
    [COLOR=#804040][b]end do[/b][/color]
    [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Result for given macheps = '[/color], macheps
    [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'t = '[/color], t,[COLOR=#ff00ff]':         sum(1/t) = '[/color], s
  [COLOR=#a020f0]end subroutine[/color] harmonic_series
[COLOR=#a020f0]end module[/color] my_numeric

[COLOR=#a020f0]program[/color] sum_of_series
  [COLOR=#a020f0]use[/color] my_numeric

[COLOR=#2e8b57][b]  real[/b][/color] :: macheps, s, s_old, t

  [COLOR=#0000ff]! compute Machine Epsilon[/color]
  macheps [COLOR=#804040][b]=[/b][/color] machine_epsilon([COLOR=#ff00ff]0.5[/color])
  [COLOR=#804040][b]write[/b][/color] ([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Machine Epsilon ='[/color], macheps

  [COLOR=#0000ff]! compute largest possible t for sum(1/t) aka harmonic series [/color]
  [COLOR=#a020f0]call[/color] harmonic_series(macheps)
[COLOR=#a020f0]end program[/color] sum_of_series

Compiling and running

$ g95 macheps.f90 -o macheps

$ macheps
 Machine Epsilon = 5.421011E-20
 Result for given macheps =  5.421011E-20
 t =  2097152. :         sum(1/t) =  15.403683

If you want to see how it works uncomment the write stataments in the function and the subroutine.

 

[mikrom]

Sorry, in my previous post is a typo - it should be:
...
machine epsilon, which is a number which fulfills the condition: macheps + 1 = 1 (analogical to 0+1 = 1)
...