#78: Problem with derivatives of power function at point zero ---------------------------+------------------------------------------------ Reporter: sebastien | Owner: sebastien Type: bug | Status: new Priority: major | Milestone: 4.2 Component: Preprocessor | Version: 4.1.0 Resolution: | Keywords: ---------------------------+------------------------------------------------ Changes (by sebastien):

* milestone: => 4.2

#78: Problem with derivatives of power function at point zero ---------------------------+------------------------------------------------ Reporter: sebastien | Owner: houtanb Type: bug | Status: accepted Priority: major | Milestone: 4.2 Component: Preprocessor | Version: 4.1.0 Resolution: | Keywords: ---------------------------+------------------------------------------------ Changes (by houtanb):

* status: assigned => accepted

#78: Problem with derivatives of power function at point zero ---------------------------+------------------------------------------------ Reporter: sebastien | Owner: houtanb Type: bug | Status: reopened Priority: major | Milestone: 4.2 Component: Preprocessor | Version: 4.1.0 Resolution: | Keywords: ---------------------------+------------------------------------------------ Changes (by sebastien):

* status: closed => reopened * resolution: fixed =>

Old description:

There is currently a problem with the derivatives of the power function f(x)=x^p^ at the point x=0.

The power function f(x)=x^p^ well defined for x>0 and for any p, since f(x)=e^p*ln(x)^.

The domain of definition of this function can be extended to x=0 with the following rules when p>0:

• f(0)=0
• f is differentiable k times at x=0, where k is the integer part of p,

and f'^k^=0

And when p is an strictly positive integer, the function is differentiable at any order, and its derivative at zero is zero.

When p is a constant, the preprocessor computes the right value of the derivatives at x=0, using simplifications rules.

But when p is a parameter, the preprocessor will not give the right value for integer values of p.

Currently the derivation rule is applied by the preprocessor, the k-th derivative of f(x) is equal to: p*(p-1)*...*(p-k+1)*x^p-k^

So when p is an integer, the (p+1)-th derivative computed at x=0 will give a NaN (0 times infinity), while it ought to be zero.

This problem typically arises in models with adjustment costs where the curvature (the exponent) of the cost is a integer parameter, and the cost is zero at steady state.

The solution is probably to create, in the preprocessor, a new operator called "derivative of power function at order k". This operator would output a code like that in the M-file:

{{{ dxp=deriv(x,p,k) % The k-th derivative of x^p

if (p > 0) && (abs(p - round(p)) < 1e-12 && k >= p dxp = 0; else dxp = x^(p-k); for i=1:k-1 dxp = dxp*p; p = p-1; end end

}}}

New description:

There is currently a problem with the derivatives of the power function f(x)=x^p^ at the point x=0.

The power function f(x)=x^p^ well defined for x>0 and for any p, since f(x)=e^p*ln(x)^.

The domain of definition of this function can be extended to x=0 with the following rules when p>0: * f(0)=0 * f is differentiable k times at x=0, where k is the integer part of p, and f'^k^=0

And when p is an strictly positive integer, the function is differentiable at any order, and its derivative at zero is zero.

When p is a constant, the preprocessor computes the right value of the derivatives at x=0, using simplifications rules.

But when p is a parameter, the preprocessor will not give the right value for integer values of p.

Currently the derivation rule is applied by the preprocessor, the k-th derivative of f(x) is equal to: p*(p-1)*...*(p-k+1)*x^p-k^

So when p is an integer, the (p+1)-th derivative computed at x=0 will give a NaN (0 times infinity), while it ought to be zero.

This problem typically arises in models with adjustment costs where the curvature (the exponent) of the cost is a integer parameter, and the cost is zero at steady state.

The solution is probably to create, in the preprocessor, a new operator called "derivative of power function at order k". This operator would output a code like that in the M-file:

{{{ dxp=deriv(x,p,k) % The k-th derivative of x^p

if abs(x) < 1e-12 && (p > 0) && (abs(p - round(p)) < 1e-12 && k >= p dxp = 0; else dxp = x^(p-k); for i=1:k-1 dxp = dxp*p; p = p-1; end end

}}}

--

Comment:

Fixed in 35450a292b99a7231219e44e0a1c7acf62858853 for M and C code.

Still needs to be implemented in bytecode.

#78: Problem with derivatives of power function at point zero ---------------------------+------------------------------------------------ Reporter: sebastien | Owner: ferhat Type: bug | Status: closed Priority: major | Milestone: 4.2 Component: Preprocessor | Version: 4.1.0 Resolution: fixed | Keywords: ---------------------------+------------------------------------------------ Changes (by sebastien):

* status: assigned => closed * resolution: => fixed

Comment:

Implemented in bytecode in 129d404eb86b845d26621f967705d55cf6fa8b40