[GiNaC-list] Issue with derivative of "abs"

[Pierangelo Masarati]

Dear GiNaC users,

we have been using GiNaC for years now, to provide our multibody 
dynamics solver, MBDyn <http://www.mbdyn.org/>, some symbolic 
differentiation capabilities (perhaps not the most appropriate use, but 
some users seem to like it).

I'll make the story as short as possible. In one of the most interesting 
applications, the user needs to:

- define the variables that represent strains (and strain rates),

- the expression for the configuration-dependent generalized forces

in short, a constitutive law. The expression's partial derivatives with 
respect to the strains and strain rates are then automatically computed 
by GiNaC, and, when needed, the numerical values are substituted and the 
expressions are evaluated.

The related code snippet looks like

>         GiNaC::symbol gEps("v");  // parameter symbols
>         GiNaC::ex gExpr;  // expressions
>         GiNaC::ex gExprDEps;  // derivatives
>
>         GiNaC::lst l;
>
>         l.append(gEps);
>         gExpr = GiNaC::ex("v*abs(v)", l);
>
>         gExprDEps = gExpr.diff(gEps);

One user recently tried something like this:

     var = v

     expr = v*abs(v)    # a turbulent viscous law

the result was

     dexpr/dvar = abs(v)+1/2*v*(v+conjugate(v))*abs(v)^(-1)

which fails when v = 0 because of the division by zero.  I tried with 
"realsymbol" (since our values can only be real); things, however, 
didn't change much:

     dexpr/dvar = v^2*abs(v)^(-1)+abs(v)

I realized that the issue is in the derivative of the abs() function; 
indeed, the proposed formula is rather general (despite failing when 
numerically evaluated with 0 as the argument).  When the argument is 
real, a more robust representation would be "d(abs(v))/dv = sign(v)".  
In the case at hand, "d(v*abs(v))/dv = 2*abs(v)" or, as a second choice, 
"d(v*abs(v))/dv = abs(v) + v*sign(v)".

First, I considered the possibility to implement the "sign()" function 
in GiNaC, but I am not familiar enough with its internals. If it makes 
sense, I would definitely recommend its introduction.

However, I noticed that the "step()" function is implemented (and 
step(0) correctly evaluates to 1/2); the "sign()" function can then be 
obtained as "sign(x) = 2*step(x) - 1" when "x" is real. To this end, I 
modified GiNaC's code as reported here:

> diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
> index d68afbb5..ec74a435 100644
> --- a/ginac/inifcns.cpp
> +++ b/ginac/inifcns.cpp
> @@ -321,6 +321,8 @@ static ex abs_expand(const ex & arg, unsigned options)
>  static ex abs_expl_derivative(const ex & arg, const symbol & s)
>  {
>         ex diff_arg = arg.diff(s);
> +       if (arg.info(info_flags::real))
> +               return diff_arg*(2*step(arg) - 1);
>         return 
> (diff_arg*arg.conjugate()+arg*diff_arg.conjugate())/2/abs(arg);
>  }
which indeed gives the expected result:

     dexpr/dvar = v*(-1+2*step(v))+abs(v)

I wonder whether my fix can be considered correct and, in case, whether 
it could be included in the distribution.

Sincerely, p.

-- 
Pierangelo Masarati
Professore Ordinario di Costruzioni e Strutture Aerospaziali
Dipartimento di Scienze e Tecnologie Aerospaziali
Politecnico di Milano