Re: Strange behavior in algebric computations

Hi Mary,

Here are the results of your codes:

ncl 0> su = 0.
ncl 1> dh = 1271
ncl 2> v1 = 0.699997
ncl 3> v2 = -4.8
ncl 4> su1 = su + dh * (v1 + v2) / 2
ncl 5> su2 = su + 0.5 * dh * (v1 + v2)
ncl 6> print("su1 = " + su1)
(0) su1 = -2605.55
ncl 7> print("su2 = " + su2)
(0) su2 = -2605.55

ncl 10> su = -2605.55
ncl 11> dh = 2114
ncl 12> v1 = -4.8
ncl 13> v2 = -10.7
ncl 14> su1 = su + dh * (v1 + v2) / 2
ncl 15> su2 = su + 0.5 * dh * (v1 + v2)
ncl 16> print("su1 = " + su1)
(0) su1 = -18989.1
ncl 17> print("su2 = " + su2)
(0) su2 = -18989.1

As I mentioned in the previous email, with both formulae the result is
the same, when computed by hand, but in the loop it goes to a missing
value.

And here are the outputs from printVarSummary for all variables used
in the formula:

Variable: dh
Type: float
Total Size: 4 bytes
            1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
Number Of Attributes: 1
  _FillValue : -32767

Variable: sdh
Type: float
Total Size: 4 bytes
            1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
Number Of Attributes: 1
  _FillValue : -32767

Variable: sv
Type: float
Total Size: 4 bytes
            1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
Number Of Attributes: 1
  _FillValue : -32767

Variable: v (subsection)
Type: float
Total Size: 4 bytes
            1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
Number Of Attributes: 8
  lon : 302.5
  lat : -35
  nivel : 250
  tempo : 1
  padrao : 1
  _FillValue : -32767
  long_name : vento meridional
  units : m/s

Variable: hgt (subsection)
Type: float
Total Size: 4 bytes
            1 values
Number of Dimensions: 1
Dimensions and sizes: [1]
Coordinates:
Number Of Attributes: 8
  lon : 302.5
  lat : -35
  nivel : 300
  tempo : 1
  padrao : 1
  _FillValue : -32767
  long_name : altura geopotencial
  units : m

And here is the loop

sdh = 0.
su = 0.
sv = 0.
dh = 0.
do k = k10000, k3000
            dh = hgt(c,t,k-1,j,i) - hgt(c,t,k,j,i)
            sdh = sdh + dh
            su = su + 0.5 * dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) )
            sv = sv + 0.5 * dh * ( v(c,t,k-1,j,i) + v(c,t,k,j,i) )
end do

And all these lines are inside other loops for c, t, i, and j. I
searched for missing values in the data (hgt, u, and v), founding
nothing. It is very strange. I thought it was due a error
accumulations, but the variables involved are equaled to zero before k
loop.

Best regards,

Mateus

2009/10/1 Mary Haley <haley_at_ucar.edu>:
> Hi Mateus,
>
> I was unable to reproduce your problem. Can you try this code and let me
> know what you get:
>
>  su = 0.
>  dh = 1271
>  v1 = 0.699997
>  v2 = -4.8
>  su1  = su +  dh * (v1 + v2) / 2
>  su2  = su +  0.5 * dh * (v1 + v2)
>  print("su1 = " + su1)
>  print("su2 = " + su2)
>
>  su = -2605.55
>  dh = 2114
>  v1 = -4.8
>  v2 = -10.7
>  su1  = su +  dh * (v1 + v2) / 2
>  su2  = su +  0.5 * dh * (v1 + v2)
>  print("su1 = " + su1)
>  print("su2 = " + su2)
>
> My first thought is that this has something to do with the types of your
> variables, and the order of which the arithmetic is being done
> on them.
>
> Can you do a "printVarSummary" on all your variables and send me the output
> so I have the correct types?
>
> --Mary
>
>
> On Oct 1, 2009, at 8:13 AM, Mateus Teixeira wrote:
>
>> Dear NCL users,
>>
>> I'm doing some algebric computations that are having strange behavior
>> depending on way the formula is written. When written this way
>>
>> su  = su +  dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) ) / 2
>>
>> I'm getting a missing value (the formula is inside a do loop), as shown
>> below:
>>
>> (0)     sv(before)=0; dh=1271; v(k-1)=0.699997; v(k)=-4.8
>> (0)     sv(after)=-2605.55
>> (0)     sv(before)=-2605.55; dh=2114; v(k-1)=-4.8; v(k)=-10.7
>> (0)     sv(after)=-32767           <= MISSING VALUE
>>
>> But when the formula is written as
>>
>> su  = su +  0.5 * dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) )
>>
>> I haven't a missingl value, as shown below for the same stop point showed
>> above
>>
>> (0)     sv(before)=0; dh=1271; v(k-1)=0.699997; v(k)=-4.8
>> (0)     sv(after)=-2605.55
>> (0)     sv(before)=-2605.55; dh=2114; v(k-1)=-4.8; v(k)=-10.7
>> (0)     sv(after)=-18989.1     <=CORRECT VALUE
>>
>> I made these same computations by hand in NCL with both formulae and
>> get the correct value given by the second way of the formula.
>>
>> Can anyone imagine what is causing this? I thought that the float
>> representation of the 1/2 could be the source of the error, but the
>> difference between missing value and the correct value is not little.
>>
>> Best regards,
>>
>> --
>> Mateus da Silva Teixeira
>> Meteorologista
>> Instituto de Pesquisas Meteorológicas - IPMet
>> Universidade Estadual Paulista - UNESP
>> Av. Luis Edmundo Carrijo Coube, 14-01 - Cx. Postal 281 - CEP 17033-360
>> Bauru - SP - Brasil
>> Fone: +55 14 3103-6030 / fax: 3203-3649
>>
>> Registered Linux User #466740 (http://counter.li.org/)
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>

-- 
Mateus da Silva Teixeira
Meteorologista
Instituto de Pesquisas Meteorológicas - IPMet
Universidade Estadual Paulista - UNESP
Av. Luis Edmundo Carrijo Coube, 14-01 - Cx. Postal 281 - CEP 17033-360
Bauru - SP - Brasil
Fone: +55 14 3103-6030 / fax: 3203-3649
Registered Linux User #466740 (http://counter.li.org/)