Unit Vector defined by two near points Help 1

[Kostas85]

Hello all,
I am trying to make a unit vector subroutine for two near points. Given the two points a1, a2 (two triplets of number) I have to define the e=(a2-a1)/distance(a1,a2) {note that a2-a1 is in vector notation}. My problem arises when I have to evaluate the above as a1 tends to a2 and more particular when lets say one has a1=(1,2,1000.00001) and a2=(1,2,1000). Then I get a NaN vector as a result. I am using single precision and I would prefer to keep it that way. Could you please propose a solution so that the computer difference result 1000.00001-1000 is actualy a 0.00001 and not a zero.
Thanks

 

[mikrom]

Hi Kostas85,

I looked at your problem. Here is an example:
unit_vector.f95

[COLOR=#a020f0]program[/color] unit_vector
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
  [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]dimension[/b][/color]([COLOR=#ff00ff]3[/color]) :: u, v, z, e

  u [COLOR=#804040][b]=[/b][/color] ([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]1.d0[/color], [COLOR=#ff00ff]2.d0[/color], [COLOR=#ff00ff]1000.00001d0[/color][COLOR=#804040][b]/[/b][/color])
  v [COLOR=#804040][b]=[/b][/color] ([COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]1.d0[/color], [COLOR=#ff00ff]2.d0[/color], [COLOR=#ff00ff]1000.d0[/color][COLOR=#804040][b]/[/b][/color])
  z [COLOR=#804040][b]=[/b][/color] u [COLOR=#804040][b]-[/b][/color] v
  e [COLOR=#804040][b]=[/b][/color] z[COLOR=#804040][b]/[/b][/color]norm(z)

  [COLOR=#0000ff]!print *, 'u = (', u, ')'[/color]
  [COLOR=#0000ff]!print *, 'v = (', v, ')'[/color]
  [COLOR=#0000ff]!print *, 'z = u - v = (', z, ')'[/color]
  [COLOR=#0000ff]!print *, 'norm(z)=', norm(z)[/color]
  [COLOR=#0000ff]!print *, 'e = (', e, ')'[/color]

  [COLOR=#804040][b]print[/b][/color] [COLOR=#ff00ff]10[/color], [COLOR=#ff00ff]'u = ('[/color], u([COLOR=#ff00ff]1[/color]), [COLOR=#ff00ff]', '[/color], u([COLOR=#ff00ff]2[/color]), [COLOR=#ff00ff]', '[/color], u([COLOR=#ff00ff]3[/color]), [COLOR=#ff00ff]')'[/color]
  [COLOR=#804040][b]print[/b][/color] [COLOR=#ff00ff]10[/color], [COLOR=#ff00ff]'v = ('[/color], v([COLOR=#ff00ff]1[/color]), [COLOR=#ff00ff]', '[/color], v([COLOR=#ff00ff]2[/color]), [COLOR=#ff00ff]', '[/color], v([COLOR=#ff00ff]3[/color]), [COLOR=#ff00ff]')'[/color]
  [COLOR=#804040][b]print[/b][/color] [COLOR=#ff00ff]10[/color], [COLOR=#ff00ff]'z = u - v = ('[/color], z([COLOR=#ff00ff]1[/color]), [COLOR=#ff00ff]', '[/color], z([COLOR=#ff00ff]2[/color]), [COLOR=#ff00ff]', '[/color], z([COLOR=#ff00ff]3[/color]), [COLOR=#ff00ff]')'[/color]
  [COLOR=#804040][b]print[/b][/color] [COLOR=#ff00ff]20[/color], [COLOR=#ff00ff]'norm(z) =  '[/color], norm(z)
  [COLOR=#804040][b]print[/b][/color] [COLOR=#ff00ff]10[/color], [COLOR=#ff00ff]'e = ('[/color], e([COLOR=#ff00ff]1[/color]), [COLOR=#ff00ff]', '[/color], e([COLOR=#ff00ff]2[/color]), [COLOR=#ff00ff]', '[/color], e([COLOR=#ff00ff]3[/color]), [COLOR=#ff00ff]')'[/color]

  [COLOR=#6a5acd]10[/color] [COLOR=#804040][b]format[/b][/color](a15 [COLOR=#008080]f10.5[/color] a2 [COLOR=#008080]f10.5[/color] a2 [COLOR=#008080]f10.5[/color] a)
  [COLOR=#6a5acd]20[/color] [COLOR=#804040][b]format[/b][/color](a15 [COLOR=#008080]f10.5[/color])

[COLOR=#a020f0]contains[/color]

[COLOR=#a020f0]function[/color] norm(vector)
  [COLOR=#0000ff]! Find the norm of vector[/color]
  [COLOR=#2e8b57][b]double precision[/b][/color] :: norm  
  [COLOR=#2e8b57][b]double precision[/b][/color], [COLOR=#2e8b57][b]intent[/b][/color] ([COLOR=#2e8b57][b]in[/b][/color]), [COLOR=#2e8b57][b]dimension[/b][/color] (:) :: vector
  [COLOR=#2e8b57][b]integer[/b][/color] :: i
  norm[COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]0[/color]
  [COLOR=#804040][b]do[/b][/color] i [COLOR=#804040][b]=[/b][/color][COLOR=#ff00ff]1[/color], [COLOR=#008080]size[/color](vector)
    norm[COLOR=#804040][b]=[/b][/color]norm[COLOR=#804040][b]+[/b][/color]vector(i)[COLOR=#804040][b]**[/b][/color][COLOR=#ff00ff]2[/color]
  [COLOR=#804040][b]end do[/b][/color]
  norm[COLOR=#804040][b]=[/b][/color] [COLOR=#008080]sqrt[/color](norm)
[COLOR=#a020f0]end function[/color] norm

[COLOR=#a020f0]end program[/color] unit_vector

I used double precision and it works fine

$ g95 unit_vector.f95 -o unit_vector

$ unit_vector
          u = (   1.00000,    2.00000, 1000.00001)
          v = (   1.00000,    2.00000, 1000.00000)
  z = u - v = (   0.00000,    0.00000,    0.00001)
    norm(z) =     0.00001
          e = (   0.00000,    0.00000,    1.00000)

When I replace in the above code double precision with single precision real I get this result

$ g95 unit_vector_real.f95 -o unit_vector_real

$ unit_vector_real
          u = (   1.00000,    2.00000, 1000.00000)
          v = (   1.00000,    2.00000, 1000.00000)
  z = u - v = (   0.00000,    0.00000,    0.00000)
    norm(z) =     0.00000
          e = (       NaN,        NaN,        NaN)

It seems that 1000.00001 could not be accurate represented using single precision float data type, instead it's represented as 1000.00000.
This causes that the computed norm will be 0 and then the division by zero delivers NaN (Not a Number)

To avoid this problem, you need to use double precision.

 

[ArkM]

Some addition:
32-bit floats keep ~7 valued decimal digits only. 1000.00001 has 9 digits...