axiom-developer
[Top][All Lists]

## [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.

----------- 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))) );

--Runs out of memory.
'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 ));
--   '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 ));
```
```

----------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.

```