- GiNaC-list - ginac.de

Re: [GiNaC-list] Issue with derivative of "abs"
by Vladimir V. Kisil 16 Jun '20

16 Jun '20

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

1 0

Release 1.7.9
by Richard B. Kreckel 11 Apr '20

11 Apr '20

Hi, GiNaC 1.7.9 has just been released. The changes are: * Fix unarchiving of overloaded functions. * Fix MinGW64 build. Happy Easter! -richy. -- Richard B. Kreckel <https://in.terlu.de/~kreckel/>

1 0

archive/unarchive on functions with a variable number of arguments
by Feng Feng 08 Apr '20

08 Apr '20

Hi all, I came across a problem about the archive and unarchive on functions with a variable number of arguments. for example, the following c++ code produces the output: zeta(10,x) :> 1/93555*Pi^10. it seems that the unarchive just return zeta(10). Best regards! Feng //--- symbol x("x"); ex zx = zeta(10,x); archive ar_out; ar_out.archive_ex(zx, "zx"); ofstream out("tmp.gar"); out << ar_out; out.close(); lst syms = {x}; archive ar_in; ifstream in("tmp.gar"); in >> ar_in; in.close(); cout << zx << " :> " << ar_in.unarchive_ex(syms,"zx") << endl; //---

3 4

Re: [GiNaC-list] series of x^n around x=0
by Feng Feng 23 Dec '19

23 Dec '19

Of course, the issue here is mathematics, not GiNaC. I totally agree with that. :) And thanks a lot for pointing out the method to find the leading power. Since there is also log(x)-terms in my expression, and leading power seems not be able to capture these log terms. Currently, I am trying to walk around by finding the all pattern: x^n, since n is rational number, and multiply n by lcm of all those denominators, and transform all rational powers to integer powers … BTW, when the input contains x^n (n as an integer) and log(x), it seems that series works fine around x=0. Thanks again! Best regards! Feng 原始邮件发件人: Vladimir V. Kisil<kisilv(a)maths.leeds.ac.uk> 收件人: Feng Feng<F.Feng(a)outlook.com> 抄送: GiNaC discussion list<ginac-list(a)ginac.de>; Vladimir V. Kisil<V.Kisil(a)leeds.ac.uk> 发送时间: 2019年12月22日(周日) 06:17 主题: Re: [GiNaC-list] series of x^n around x=0 >>>>> On Sat, 21 Dec 2019 11:14:58 +0000, Feng Feng <F.Feng(a)outlook.com<mailto:F.Feng@outlook.com>> said: FF> Thanks very much for the reply.! And yes, it is not FF> well-defined in mathematical sense to taylor expand x^(3/2) FF> around x=0. What I exactly want is to get the asyptotic FF> behaviour of a function f(x) around x=0. For example, f(x) = FF> x^(3/2) (1 + x + x^2 + …), the asyptotic behaviour at leading FF> order (LO) and next-to-leading order (NLO) are x^(3/2) and FF> x^(5/3), respectively. FF> So I wonder there is a way to get the asyptotic behaviour of a FF> function f(x), for example, Input: x^(3/2) * (1+x+x^2+x^3), and FF> request the terms at LO and NLO Output: x^(3/2) + x^(5/3) I am not sure that is NLO, but the leading term can be obtained as follows: ex e=pow(x, numeric(3,2))*(1+x+x^2+x^3); cout << (log(e).diff(x)*x).expand().subs(x==0); // -> 3/2 So now e can be multiplied by x^(2/3) to remove the singularity. See for explanations: https://en.wikipedia.org/wiki/Logarithmic_derivative So, the issue here is mathematics, not GiNaC. Once you know the right method, it shall be easier to implement in the code. -- 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 (Colab): https://github.com/vvkisil/MoebInv-notebooks Jupyter (CodeOcean): https://codeocean.com/capsule/7952650/tree

1 0

Re: [GiNaC-list] series of x^n around x=0
by Feng Feng 21 Dec '19

21 Dec '19

Thanks very much for the reply.! And yes, it is not well-defined in mathematical sense to taylor expand x^(3/2) around x=0. What I exactly want is to get the asyptotic behaviour of a function f(x) around x=0. For example, f(x) = x^(3/2) (1 + x + x^2 + …), the asyptotic behaviour at leading order (LO) and next-to-leading order (NLO) are x^(3/2) and x^(5/3), respectively. So I wonder there is a way to get the asyptotic behaviour of a function f(x), for example, Input: x^(3/2) * (1+x+x^2+x^3), and request the terms at LO and NLO Output: x^(3/2) + x^(5/3) Thanks again! Best regards! Feng

2 1

series of x^n around x=0
by Feng Feng 19 Dec '19

19 Dec '19

Hi all, I got a question while using series on x^n with n a non-integer. here is an example: ------ symbol x("x"); ex expr = pow(x, ex(3)/2); // 3/2 may be a general non-integer cout << expr.series(x,2) << endl; ------ The output is: Order(x^2) so is there a way to handle this? Best regards! Feng

2 1

Undefine the predefined numeric consant 'I'?
by Warren Weckesser 03 Nov '19

03 Nov '19

Is there an easy way to undefine the predefined numeric value 'I'? If I define a symbol 'I', use it in an expression, and then output that expression to cout after modifying stream with csrc, 'I' is converted to 'std::complex<double>(0.0,-1.0)'. I'd like the output to be just 'I'. Warren

3 4

Question: "expand" power in ginac
by Clobert von Roum 16 Oct '19

16 Oct '19

3 3

Release 1.7.8
by Richard B. Kreckel 07 Oct '19

07 Oct '19

Hi, GiNaC 1.7.8 has just been released. The changes are: * Fix pseries::evalf(), broken since 1.7.0. * Fix a corner-case bug in H_evalf(). Happy hacking! -richy. -- Richard B. Kreckel <https://in.terlu.de/~kreckel/>

1 0

Release 1.7.7
by Richard B. Kreckel 22 Sep '19

22 Sep '19

Hi, GiNaC 1.7.7 has just been released. It fixes some bugs (see NEWS) and depends on Python3 for build (this may be important for distros). Happy hacking! -richy. -- Richard B. Kreckel <https://in.terlu.de/~kreckel/>

1 0