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