axiom-developer
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Axiom-developer] Higher order derivatives.


From: A.B.
Subject: [Axiom-developer] Higher order derivatives.
Date: Mon, 01 Apr 2013 19:15:22 +0400
User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:17.0) Gecko/20130311 Thunderbird/17.0.4

Hello.

I have learnt about axiom recently and now I'm trying to apply it to some of my current problems (numerical computations in physics).

I need to evaluate higher order derivatives of Lewis integral ( http://link.aps.org/doi/10.1103/PhysRev.102.537 , Appendix A, eq. 9) for further fortran export.

Using a straightforward approach (see -----v1 below), I run out of memory already on the second derivative.

Following an example from the 'Derivatives' section of the axiom book (sec.1.11), I can evaluate the derivative I need in a general form. (see ----v2). But when I substitute actual expressions for operators LiBt, LiDt and LiAg, I run out of memory again.

Another way (see ---v3) could be to replace derivatives of operators in the general form by some dummy variables, then actually define them and export to fortran. This way looks promising but I do not know how to automate it, so it becomes impractical for derivatives of 2nd order and higher.

So, my question is: what is the best way to evaluate, for example, 5th derivative of 'Li' function from examples below and
export this result to fortran?

P.S: sorry for my english.

--------------- v1

)clear all

LiBt :=_
  mu * ((qx-px)^2 + (qy-py)^2 + (qz-pz)^2 + (a+b)^2 )_
  + b * ( mu^2 + qx^2 + qy^2 + qz^2 + a^2 )_
  + a * ( mu^2 + px^2 + py^2 + pz^2 + b^2 );

LiAg :=_
  ( (qx-px)^2 + (qy-py)^2 + (qz-pz)^2 + (a+b)^2 )_
  * ( qx^2 + qy^2 + qz^2 + (mu+a)^2 )_
  * ( px^2 + py^2 + pz^2 + (mu+b)^2 );

LiD := sqrt( LiBt^2 - LiAg );

Li := %pi^2 / LiD * log( (LiBt + LiD) / (LiBt - LiD) );

--Runs out of memory during computation.
dLidMudAdB := D(Li,[mu,a,b]);

----------- v2

)clear all

LiDt := operator 'LiDt
LiBt := operator 'LiBt
LiAg := operator 'LiAg

Li := (%pi^2 / LiDt(LiBt(mu,a,b),LiAg(mu,a,b)))_
    * log(_
     ( LiBt(mu,a,b) + LiDt( LiBt(mu,a,b), LiAg(mu,a,b)))_
     / ( LiBt(mu,a,b) - LiDt( LiBt(mu,a,b), LiAg(mu,a,b))) );

dLidMudAdBGeneral := D(Li,[mu,a,b]);

--Runs out of memory.
dLidMudAdBsubstBt := eval(dLidMudAdBGeneral,_
  'LiBt, (x) +-> x.1 * ( (qx-px)^2 + (qy-py)^2 + (qz-pz)^2 + (x.2+x.3)^2 )_
  + x.3 * ( x.1^2 + qx^2 + qy^2 + qz^2 + x.2^2 )_
  + x.2 * ( x.1^2 + px^2 + py^2 + pz^2 + x.3^2 ));
-- dLidMudAdBsubstAg := eval(dLidMudAdBsubstBt,_
--   'LiAg, (x) +-> ( (qx-px)^2 + (qy-py)^2 + (qz-pz)^2 + (x.2+x.3)^2 )_
--   * ( qx^2 + qy^2 + qz^2 + (x.1+x.2)^2 )_
--   * ( px^2 + py^2 + pz^2 + (x.1+x.3)^2 ));
-- dLidMudAdBsubstDt := eval(dLidMudAdBsubstAg, 'LiDt, (x) +-> sqrt(x.1^2 - x.2) );


----------v3

)clear all

LiDt := operator 'LiDt
LiBt := operator 'LiBt
LiAg := operator 'LiAg

Li := (%pi^2 / LiDt(LiBt(mu,a,b),LiAg(mu,a,b)))_
    * log(_
     ( LiBt(mu,a,b) + LiDt( LiBt(mu,a,b), LiAg(mu,a,b)))_
     / ( LiBt(mu,a,b) - LiDt( LiBt(mu,a,b), LiAg(mu,a,b))) )

dLidMuGeneral := D(Li,mu)

d1DRules := rule
    D(LiDt(beta,alphagamma),[beta]) == dDdB
    D(LiDt(beta,alphagamma),[alphagamma]) == dDdAg
d0DRules := rule
    LiDt(beta,alphagamma) == LiD
d1BRules := rule
    D(LiBt(mu,a,b),[mu]) == dBdMu
d0BRules := rule
    LiBt(mu,a,b) == LiB
d1AgRules := rule
    D(LiAg(mu,a,b),[mu]) == dAgdMu
d0AgRules := rule
    LiAg(mu,a,b) == LiA

dLidMuSubst := d0AgRules(d0BRules(d0DRules(d1AgRules(d1BRules(d1DRules(dLidMuGeneral))))))
-- Then export dLidMuSubst to fortran.
-- Then define
--   LiD := sqrt( b^2 - ag)
--   dDdB := D(LiD,b)
--   dDdA := D(LiD,ag)
-- and similar definitions for dBdMu, dAgdMu, LiB, LiA.
-- Then somehow export them.




reply via email to

[Prev in Thread] Current Thread [Next in Thread]