a.a.letterbox wrote:
>
> I need to evaluate higher order derivatives of Lewis integral (
> 10.1103/PhysRev.102.537 , Appendix A, eq. 9) for further fortran export.
>
> With a straightforward approach (see -----v1 below), I'm running out of
> memory already on the second derivative.
>
> I was able to progress slightly further (see ----v2) adapting an example
> from the 'Derivatives' chapter of the axiom book (sec.1.11), but
> i don't know how to export results into fortran, since after evaluation of
> dLidMudAdB I need to somehow define and export all functions like
> D(bt(mu,a,b),[a,b,mu]),
> D(bt(mu,a,b),[b,mu]) and so on.
>
> All in all, my question is: what is the best way to compute, say, 5th
> derivative of 'Li' function from examples below and
> export this result to fortran?
>
> P.S: I was under impression that attachements to mail lists are not a good
> idea.
> If they are fine actually, next time i'll try to use them to include code.
Plain text is preferable.
>
> 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) );
>
> dLidMudA := D(Li,[mu,a]);
>
You have not so small expression for Li. Derivatives are
going to be much bigger. Axiom and FriCAS keeps expressions
in fully expanded form, so the size is likely to overflow
memory. Even if some derivatives can be computed, the
expression probably will be too large for efficient
numerical exaluation.
To efficienty compute derivatives one needs to reuse
subexpressions: derivatives have _a lot_ of repeated
subexpressions. Basically one need something like
"stright line program". In Fortran one wants to
store value of each subexpression in a variable
and use this variable in further computations.
Your v2 is a step in this direction. Below I
give an example of computing second derivative
(sorry, third derivative would require more
work and second already explains the idea).
In this example I store values of LiBt, LiAg
and their derivatives. Storing values and
derivatives of LiD (as a function of LiBt and
LiAg) would probably give smaller program, but
requires more effort. Basic idea is that
in the result of your v2 we need to substitute
actual values for operators. We do this using
'eval' function. We need separate substitutions
for base functions and for derivatives. In
the example I have build list of substitutions
by hand. It shoud be possible to build the
substitution list in mechanical way. However,
there is little tricky point here: we need
names which will be valid Fortran variable
names. For this one probably should first
generate names as strings and then convert
strings to symbol and then to expressions
when building substitution list.
One extra remark: I just noticed that in
current FriCAS Fortran output of exponents
is broken, I used older version to get output.
------- script below
)clear all
bt := operator 'bt
ag := operator 'ag
LiBt := bt(mu,a)
LiAg := ag(mu,a)
LiD := sqrt( LiBt^2 - LiAg );
Li := %pi^2 / LiD * log( (LiBt + LiD) / (LiBt - LiD) );
dLidMudA0 := D(Li,[mu,a]);
)clear values LiBt LiAg
dLidMudA := eval(dLidMudA0, [bt(a,b) = LiBt, ag(a,b) = LiAg, _
D(bt(mu,a), a) = LiBt_a, _
D(bt(mu,a), mu) = LiBt_mu, _
D(ag(mu,a), a) = LiAg_a, _
D(ag(mu,a), mu) = LiAg_mu, _
D(bt(mu,a), [mu, a]) = LiBt_mu_a, _
D(ag(mu,a), [mu, a]) = LiAg_mu_a]);
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 );
LiBt_a := D(LiBt, a)
LiBt_mu := D(LiBt, mu)
LiAg_a := D(LiAg, a)
LiAg_mu := D(LiAg, mu)
LiBt_mu_a := D(LiBt_mu, a)
LiAg_mu_a := D(LiAg_mu, a)
outputAsFortran("LiBta", LiBt_a)
outputAsFortran("LiBtmu", LiBt_mu)
outputAsFortran("LiAga", LiAg_a)
outputAsFortran("LiAgmu", LiAg_mu)
outputAsFortran("LiBtmua", LiBt_mu_a)
outputAsFortran("LiAgmua", LiAg_mu_a)
outputAsFortran("dLidMudA", dLidMudA)
--
Waldek Hebisch
[email protected]
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/fricas-devel?hl=en.
For more options, visit https://groups.google.com/groups/opt_out.
