Dear Pierangelo, Thanks for the email. It is nice to see another usage of GiNaC! Your patch combines the best from both worlds as was discussed: https://www.ginac.de/pipermail/ginac-devel/2014-April/002105.html So I am voting in favour of this patch to be applied! Meanwhile, it seems the faster solution for your would be to clone abs function with the alteration towards the desired derivative. It is not difficult: 1. Cut & Paste everything for definition abs from GiNaC's inifcns.h and inifcns.cpp to your code. 2. Replace the derivative function by abs_expl_derivative() from your patch. 3. Search & replace abs -> myabs in this portion of code. Then the function myabs() shall be ready to be used. Best wishes, Vladimir -- Vladimir V. Kisil http://www.maths.leeds.ac.uk/~kisilv/ Book: Geometry of Mobius Transformations http://goo.gl/EaG2Vu Software: Geometry of cycles http://moebinv.sourceforge.net/ Jupyter: https://github.com/vvkisil/MoebInv-notebooks
On Tue, 16 Jun 2020 19:24:32 +0200, Pierangelo Masarati <pierangelo.masar\ ati@polimi.it> said:
PM> Dear GiNaC users, PM> we have been using GiNaC for years now, to provide our multibody PM> dynamics solver, MBDyn <http://www.mbdyn.org/>, some symbolic PM> differentiation capabilities (perhaps not the most appropriate PM> use, but some users seem to like it). PM> I'll make the story as short as possible. In one of the most PM> interesting applications, the user needs to: PM> - define the variables that represent strains (and strain PM> rates), PM> - the expression for the configuration-dependent generalized PM> forces PM> in short, a constitutive law. The expression's partial PM> derivatives with respect to the strains and strain rates are PM> then automatically computed by GiNaC, and, when needed, the PM> numerical values are substituted and the expressions are PM> evaluated. PM> 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); PM> One user recently tried something like this: PM> var = v PM> expr = v*abs(v) # a turbulent viscous law PM> the result was PM> dexpr/dvar = abs(v)+1/2*v*(v+conjugate(v))*abs(v)^(-1) PM> which fails when v = 0 because of the division by zero. I tried PM> with "realsymbol" (since our values can only be real); things, PM> however, didn't change much: PM> dexpr/dvar = v^2*abs(v)^(-1)+abs(v) PM> I realized that the issue is in the derivative of the abs() PM> function; indeed, the proposed formula is rather general PM> (despite failing when numerically evaluated with 0 as the PM> argument). When the argument is real, a more robust PM> representation would be "d(abs(v))/dv = sign(v)". In the case PM> at hand, "d(v*abs(v))/dv = 2*abs(v)" or, as a second choice, PM> "d(v*abs(v))/dv = abs(v) + v*sign(v)". PM> First, I considered the possibility to implement the "sign()" PM> function in GiNaC, but I am not familiar enough with its PM> internals. If it makes sense, I would definitely recommend its PM> introduction. PM> However, I noticed that the "step()" function is implemented PM> (and step(0) correctly evaluates to 1/2); the "sign()" function PM> can then be obtained as "sign(x) = 2*step(x) - 1" when "x" is PM> 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); >> } PM> which indeed gives the expected result: PM> dexpr/dvar = v*(-1+2*step(v))+abs(v) PM> I wonder whether my fix can be considered correct and, in case, PM> whether it could be included in the distribution. PM> Sincerely, p. PM> -- Pierangelo Masarati Professore Ordinario di Costruzioni e PM> Strutture Aerospaziali Dipartimento di Scienze e Tecnologie PM> Aerospaziali Politecnico di Milano PM> _______________________________________________ GiNaC-list PM> mailing list GiNaC-list@ginac.de PM> https://www.ginac.de/mailman/listinfo/ginac-list