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

Tue Jun 16 21:32:56 CEST 2020

	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 at 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 at ginac.de                                        
    PM> https://www.ginac.de/mailman/listinfo/ginac-list