[gradsusr] odd behavior of pow() function

[Jeff Duda]

Eric and others,
Thanks for the insights. I totally neglected the issue that there could be
problems taking non-integer powers of negative numbers.  However, more
importantly, I discovered there was an error in my calculation that I
fixed.  Since then, the problem has been fixed for the most part.

Jeff

On Fri, Jan 18, 2013 at 7:26 PM, Eric Altshuler <ela at cola.iges.org> wrote:

> There seem to be two separate issues here. One is the pow function
> returning 'nan' instead of the missing data value. As Jennifer says, this
> is a bug. However, the other issue is that lev-e should never be negative.
> The vapor pressure, e, should be much smaller than the pressure, lev (since
> lev is in millibars, e must also be in mb). In the calculation of e, I
> suspect 'rhprs' is given as a percentage, in which case e will be a factor
> of 100 too large. Multiply by 0.0611 instead of 6.11 to fix this problem.
> (Since I'm not familiar with your other formulas, I can't vouch for their
> correctness.)
>
> Using maskout in this situation will get rid of the nan's but it won't fix
> the underlying problem. You'll get missing values where lev-e is
> (erroneously) negative, and non-missing but spurious values where lev-e is
> positive (but still wrong).
>
> Finally, your calculator's order of operations may perform the ^ before
> the negation: -(3.55595^0.2854) = -1.43629. No calculator should give a
> real result for (-3.55595)^0.2854; this is mathematically impossible. Only
> if the power is the exact reciprocal of an odd integer will the raising of
> a negative number to a (decimal) non-integer power result in a real value.
> As far as I can tell, this is only true for reciprocals of powers of 5
> (1/5, 1/25, 1/125, etc.), i.e. 0.2, 0.04, 0.008, etc. Note that the power
> 1/3 (cube root) of a negative number is real, but the decimal
> 0.333333333... doesn't represent 1/3 exactly, so some calculators will give
> an error. (On calculators that accept fractional notation, the operation
> (-8)^(1/3) does give -2, as expected, when 1/3 is given as a fraction
> instead of a decimal. This is also true for other odd roots, e.g.
> (-128)^(1/7)=-2.)
>
> Best regards,
>
>
> Eric L. Altshuler
> Research Scientist
> Center for Ocean-Land-Atmosphere Studies
>
> 4041 Powder Mill Road, Suite 302
> Calverton, MD 20705-3106
> USA
>
> E-mail: ela at cola.iges.org
> Phone: (301) 902-1257
> Fax: (301) 595-9793
>
> ------------------------------
> *From: *"Jennifer Adams" <jma at cola.iges.org>
> *To: *"GrADS Users Forum" <gradsusr at gradsusr.org>
> *Sent: *Friday, January 18, 2013 2:09:24 PM
> *Subject: *Re: [gradsusr] odd behavior of pow() function
>
>
> Interesting question! From the man page on 'pow': pow(x, y) returns a NaN
> and raises the "invalid" floating-point exception for finite x < 0 and
> finite non-integer y.
>
> Here are some GrADS calls that illustrate this functionality of 'pow':
> ga-> d pow(2,0.5)
> Result value = 1.41421
> ga-> d pow(-2,2)
> Result value = 4
> ga-> d pow(-2,-1)
> Result value = -0.5
> ga-> d pow(-2,0.9999999999)
> Result value = nan
> ga-> d pow(-2,1.1)
> Result value = nan
>
> The problem is that GrADS is returning the result of the call to 'pow'
> without checking whether it is a NaN and turning it into a proper missing
> data value. I will add this to my list of bugs to fix. Meanwhile, wrap a
> maskout around the expression 1000/(lev-e) and that should any negative
> values for the first argument to 'pow'.
> --Jennifer
>
>
>
> On Jan 18, 2013, at 1:21 PM, Jeff Duda wrote:
>
> Hello,
> I am trying to compute some thermodynamic parameters, shown below.
>
>  'define e = rhprs*6.11*exp((2.5e6/461.5)*((1/273)-(1/tmpprs)))'
>  'define r = 622*e/(lev-e)'
>  'define condtemp = 55 + 1 / (1/(tmpprs-55) - log(rhprs/100)/2840)'
>  'define thetaDL =
> tmpprs*pow(1000/(lev-e),0.2854)*pow(tmpprs/condtemp,0.0028*r)'
>  'define thetae =
> thetaDL*exp(((3.036/condtemp)-0.00178)*r*(1+0.448+r/1000))'
>  'define thetaw = 45.114 - 51.489*pow(thetae/273.15,-3.504)'
>  'define LSI = max(thetaw,lev=1000,lev=500) - max(thetaw,lev=1000,lev=900)'
>
> I get a lot of NaNs.  They are coming from this line:
>  'define thetaDL =
> tmpprs*pow(1000/(lev-e),0.2854)*pow(tmpprs/condtemp,0.0028*r)'
>
> And in fact, from this piece:
> pow(1000/(lev-e),0.2854)
>
> I know this because I tested this for a model grid column (from 1000 - 500
> mb).
> ga-> d 1000/(lev-e)
> Printing Grid -- 13 Values -- Undef = -9.99e+08
> -3.55595 -4.04332 -4.69364 -9.98891 39.3747 -5.3341 -11.5849 4.19346
> 2.07427 2.42358 2.11214 2.1287 2.24423
> ga-> d pow(1000/(lev-e),0.2854)
> Printing Grid -- 13 Values -- Undef = -9.99e+08
> -nan -nan -nan -nan 2.85282 -nan -nan 1.5055
> 1.2315 1.28743 1.23787 1.24063 1.25949
>
> For whatever reason, the pow() function doesn't like negative bases,
> because for example, I can use a calculator to compute -3.55595^0.2854 as
> -1.43629.  Why, then, does the pow() function give a -nan as a result?
>
> Jeff Duda
>
> --
> Jeff Duda
> Graduate research assistant
> University of Oklahoma School of Meteorology
> Center for Analysis and Prediction of Storms
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
>
>
> --
> Jennifer M. Adams
> IGES/COLA
> 4041 Powder Mill Road, Suite 302
> Calverton, MD 20705
> jma at cola.iges.org
>
>
>
>
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
>
> _______________________________________________
> gradsusr mailing list
> gradsusr at gradsusr.org
> http://gradsusr.org/mailman/listinfo/gradsusr
>
>


-- 
Jeff Duda
Graduate research assistant
University of Oklahoma School of Meteorology
Center for Analysis and Prediction of Storms
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20130118/48af6964/attachment-0003.html