implementation of power function in fortran
hwon...@yahoo.com
I was wondering if anyone could help me with a problem I am having
with the power function in fortran
The equation I wich to implement is of the form:
Y = A*X**n
A and n are parmaters being estimated (real*8). Y and X are the
dependent and independent variable respectively.
I keep getting invalid floating point operation as an error and
suspect that it is because X is close to 0 or negative (X is a
concentation that falls).
Does anyone have any suggestions as to how to keep X positive or an
alternative to correct this problem.
Many thanks,
HW
mecej4
> Hi
>
> I was wondering if anyone could help me with a problem I am having
> with the power function in fortran
>
> The equation I wich to implement is of the form:
>
> Y = A*X**n
>
> A and n are parmaters being estimated (real*8). Y and X are the
> dependent and independent variable respectively.
>
> I keep getting invalid floating point operation as an error and
> suspect that it is because X is close to 0 or negative (X is a
> concentation that falls).
If it is true that the concentration can never become negative, you can use
if(X.gt.0.0)then
Y=A*X**n
else
Y=0
endif
> Does anyone have any suggestions as to how to keep X positive or an
> alternative to correct this problem.
>
> Many thanks,
>
> HW
HTH
-- mecej4
Arjan
>
> > A and n are parmaters being estimated (real*8). Y and X are the
> > dependent and independent variable respectively.
Take the LOG of both sides:
LOG(Y) = LOG(A) + n*LOG(X)
and use exact relations to fit a straight line!
A.
glen herrmannsfeldt
>> > Y = A*X**n
>> > A and n are parmaters being estimated (real*8). ?Y and X are the
>> > dependent and independent variable respectively.
That is a common problem with non-linear optimization
algorithms. There is much written about it, but first you
have to make sure that you don't supply negative (REAL) values
for X.
> Take the LOG of both sides:
> LOG(Y) = LOG(A) + n*LOG(X)
> and use exact relations to fit a straight line!
The OP didn't give all the details, but in some cases that works.
Though if you do, you should change (or add) the weighting function
as it changes the relative weights of the different points.
-- glen
hwon...@yahoo.com
wrote:
Thank you all.
Keep X positive using the "if then" seems to have worked.
Ron Shepard
<894080c0-3877-4719...@y6g2000prf.googlegroups.com>,
hwon...@yahoo.com wrote:
> Thank you all.
>
> Keep X positive using the "if then" seems to have worked.
Other approaches for this are to replace "X" by max(X,0.0) or abs(X), or
to fit to the function
Y = A * (X * X)**n
rather than your original expression. This latter approach is commonly
used in graphics programs where the expression is interpreted, not so
much in a compiled language like fortran. The "n" that you get this way
is half the value of the original expression, so you have to account for
that if you use "n" elsewhere.
Just as a matter of style, I would suggest using a variable like "B"
rather than "n" for a floating point exponent.
$.02 -Ron Shepard
mecej4
> In article
> <894080c0-3877-4719...@y6g2000prf.googlegroups.com>,
> hwon...@yahoo.com wrote:
>
>> Thank you all.
>>
>> Keep X positive using the "if then" seems to have worked.
>
> Other approaches for this are to replace "X" by max(X,0.0) or abs(X), or
> to fit to the function
These occurred to me, as well. But, trying to raise zero to a power would raise
exceptions. In my opinion, using abs() here is a no-no.
Could have used max(X,a small positive number), instead. Better, still, would be
to model the physics correctly.
-- mecej4
robin
news:5daff478-6836-4481...@y6g2000prf.googlegroups.com...
>> > Y = A*X**n
>
>> > A and n are parmaters being estimated (real*8). Y and X are the
>> > dependent and independent variable respectively.
>Take the LOG of both sides:
>LOG(Y) = LOG(A) + n*LOG(X)
This makes the problem worse.
If X is negative, LOG is invalid.
If X is zero, LOG is invalid.
A simple test suffices.
fj
> "Arjan" <arjan.van.d...@rivm.nl> wrote in message
Not always : in fact, using logarithms is useful when yo do that on
the theoretical part of your problem. This is a variable change. So
the function LOG is never used after that :
The equation becomes simply :
YY = AA + n * XX
and you can get the original X using :
X=EXP(XX)
I uses sometimes this transformation when solving chemistry problems
where the variable X must cover many order of magnitudes (for instance
from 1.E-20 to 1.). Because of rounding errors, it is easy to observe
that the variable XX is often more convenient (from -46. to 0.), and
the equations above cannot raise severe exceptions (except underflow
with EXP).
Ron Shepard
<7872e0a0-731c-41a5...@l5g2000vbc.googlegroups.com>,
fj <franco...@irsn.fr> wrote:
I guess this might be an obvious consequence of the change of
coordinates, but if you do least squares fitting of data to the
function
Y = a * X ** B
and to the function
YY = AA + B * XX
where the axes are transformed as YY=LOG(Y) and XX=LOG(X), then you
will get different fits because the deviations are weighted
differently.
In both cases, zero and negative values of X will cause problems,
and therefore must be accounted for. In fact, I think in fortran
that X**B is defined to be exp(B*log(X)) in the "as if" way things
like this are defined. So negative X causes problems in the two
fitting approaches for the same underlying reason, the log(X)
argument is out of range.
$.02 -Ron Shepard
Richard Maine
> In fact, I think in fortran
> that X**B is defined to be exp(B*log(X)) in the "as if" way things
> like this are defined.
No, unless I misunderstand you. It is defined simply as exponentiation,
with the elaboration that it is X raised to the B power and the
limitation that X be non-negative. It can plausibly be implemented as
exp(B*log(X)) (presumably with some special case for X being 0.0), so
perhaps that's what you meant, but that is not the definition.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
glen herrmannsfeldt
(someone wrote)
>> > >> > Y = A*X**n
>> > >> > A and n are parmaters being estimated (real*8). Y and X are the
>> > >> > dependent and independent variable respectively.
>> > >Take the LOG of both sides:
>> > >LOG(Y) = LOG(A) + n*LOG(X)
(snip)
> I guess this might be an obvious consequence of the change of
> coordinates, but if you do least squares fitting of data to the
> function
> Y = a * X ** B
> and to the function
> YY = AA + B * XX
> where the axes are transformed as YY=LOG(Y) and XX=LOG(X), then you
> will get different fits because the deviations are weighted
> differently.
That is easily fixed by doing a weighted fit. It seems rare
for anyone to do that, though. You should at least understand
the difference.
> In both cases, zero and negative values of X will cause problems,
> and therefore must be accounted for. In fact, I think in fortran
> that X**B is defined to be exp(B*log(X)) in the "as if" way things
> like this are defined. So negative X causes problems in the two
> fitting approaches for the same underlying reason, the log(X)
> argument is out of range.
There is a general problem with doing non-linear least squares
fits on functions with a restricted domain. The simple solution
is to force the variable in question back in range, but that doesn't
work so well. Better is to put a penalty with a sharply rising
but not infinitely sharp penalty function based on how far
negative (in this case) the value is. That allows the minimization
algorithm to try to correct the problem.
Any good reference on non-linear least-squares fits should
explain this.
Not knowing what the OP is actually trying to do, there is some
discussion in "Deconvolution of Images and Spectra" by Jansson.
(It seems to be unusually expensive, find one in a library.)
Jansson explains the problem of doing deconvolution on absorption
spectra, which must be between 0 and 1, with linear deconvolution
algorithms. Then he explains non-linear algorithms and how
they get around that problem.
-- glen
Ron Shepard
nos...@see.signature (Richard Maine) wrote:
> Ron Shepard <ron-s...@NOSPAM.comcast.net> wrote:
>
> > In fact, I think in fortran
> > that X**B is defined to be exp(B*log(X)) in the "as if" way things
> > like this are defined.
>
> No, unless I misunderstand you. It is defined simply as exponentiation,
> with the elaboration that it is X raised to the B power and the
> limitation that X be non-negative. It can plausibly be implemented as
> exp(B*log(X)) (presumably with some special case for X being 0.0), so
> perhaps that's what you meant, but that is not the definition.
Ok, my memory must be wrong on this. I see this definition for
complex variables (e.g. in 6.1.4 of the F77 standard), but I don't
see it anywhere defined for real variables.
$.02 -Ron Shepard
Franken Sense
I'm just always impressed with the way the fortran implementations I use
(gfortran, g95, silverfrost) get x**y (x,y real) correct to its full width
and quickly. I never forget this.
Glen brought up i to the i a while back. I think this is multi-valued, and
I couldn't think of a way to catch the other solutions.
--
Frank
If you put the two Bushs together in their over seven years of their two
presidencies, not one new job has been created. Numbers do not lie. If you
extrapolated from that, if the Bushs had run this country from its very
beginning to the current time, not one American would have ever worked.
We'd be hunter-gatherers.
~~ Al Franken, in response to the 2004 SOTU address