[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.
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Axiom-developer] Higher order derivatives.,
A.B. <=