A general template file for a Runge-Kutta function subroutine for the body
problem is detailed in this section. This involves evaluating the first
derivatives of the Hamiltonian. The subroutine below is written so that
derivatives with respect to components of position are evaluated first, followed
by derivatives of components of momentum.
In[4]:= !!rksub.mf
subroutine evalf(f,q,p)
implicit double precision (a-h,o-z)
dimension q(<*
(***** Input degrees of freedom (dimension/2) *****)
d = 3;
(***** End of required input *****)
(* Ensure that arrays are protected form N during translation. *)
SetOptions[FortranAssign,AssignToArray->{q,p}];
(* Automatically evaluate variables and n-body Hamiltonian *)
qvars = Array[q,d]; pvars = Array[p,d];
H = pvars.pvars/2 + 1/Sqrt[qvars.qvars];
d *>),p(<* d *>),f(<* 2 d *>)
<* (***** Calculate derivatives of the Hamiltonian. *****)
FortranAssign[
f,
Flatten[{Outer[D,{H},pvars],- Outer[D,{H},qvars]}]
] *>
return
end
Appropriate array dimensions are specified by Mathematica. This means that precise
control can be exercised over the generated FORTRAN subroutine, by allocating
only as much stack space as the problem dictates. This is useful for large scale
problems such as molecular dynamics simulations.
An attraction of this approach is that it may be generalised to account for any functional form of Hamiltonian. In order to do so, the Hamiltonian must be specified by the user in the template file, along with the number of degrees of freedom. Mathematica can then be used to Splice the results into a file containing FORTRAN source code for the problem.
In[5]:= Splice["rksub.mf",FormatType->OutputForm];
In[6]:= !!rksub.f
subroutine evalf(f,q,p)
implicit double precision (a-h,o-z)
dimension q(3),p(3),f(6)
f(1)=p(1)
f(2)=p(2)
f(3)=p(3)
f(4)=q(1)/sqrt((q(1)**2+q(2)**2+q(3)**2)**3)
f(5)=q(2)/sqrt((q(1)**2+q(2)**2+q(3)**2)**3)
f(6)=q(3)/sqrt((q(1)**2+q(2)**2+q(3)**2)**3)
return
end