A description of the code used to locate the minimum points of the action will be given in this subsection, including some of its specifications. The code, inventively called AVP, was written in the Fortran 90/95 programming language, consisting of subroutines, functions, and modules. The complete code will be given in Appendix A, where also a more detailed description is given of the structure of the code dependent on the programming language.
Before looking into the details of the AVP code, it should be noted that all parameters and variables are scaled according to mass, distance, and time given in the following units
`` 
'' indicating solar value. The reason for this scaling is to make
the values used in the calculation as close to unity as possible, hence reducing
the risk of roundoff errors. The resulting range of the variables, in this case the
coefficients 
in equation (3.21), is about 
.
For some of the optimizing methods this scaling is not adequate, especially in the case
of the Steepest Descent. Attempts to improve the scaling have not been successful since
there are no constraints on the value of each separate coefficient.
The computation is initialized by reading certain parameters from input files, e.g. cosmological and mass tracer parameters. The mass of each mass tracer is scaled to the mass of M31, its distance from the MW is given in Mpc, and the celestial coordinates in degrees. The celestial coordinates and distances are used to calculate each mass tracer's Cartesian coordinates by the following relations:
where d is the distance, 
is the declination and 
is the
right ascension. The mass of each mass tracer is determined by an initial guess for the
mass of M31, whereupon the celestial coordinates are adjusted to the centre of mass rest
frame.
Also read from file are the initial values of the coefficients, i.e. the starting point
for the search. As discussed earlier, this point is chosen as close as possible to
where the solution is assumed to be, which, hopefully, will reduce the computing time.
If one can assume that the system of mass tracers is cosmologically young; i.e. they are
moving towards each other for the first time, the peculiar velocities are probably small.
This indicates, by referring to equation (3.21), small values of the
coefficients. Therefore, a preferable choice for initial values should be close to zero,
both positive and negative. A common choice is actually 
. Details on
the input files are given in Appendix A.4.
Usually the coefficients are given a common initial value, but it is sometimes practical to be able to specify a specific value for each coefficient. Hence, in the AVP code, one has the choice of either giving the same value to each coefficient, or reading individual values from file. This also gives us the freedom to halt and restart the computation at any point or to use the solution from one cosmological model as initial values for a search in another.
The next step is the initialization of some parameters and variables, e.g. 
from
equation (3.36), and the trial functions (either 3.22 or
3.23). The latter ones are determined at the points in time, or rather
points in the expansion parameter, determined by the integration method. As stated in
Section 4.2, the integration method used here is the Gaussian Quadrature
method, and the number of abscissas is consistently chosen to exceed the number of trial
functions, N, by at least one. The abscissas and their respective weights are
calculated as shown in equation (4.4), where the Gauss-Legendre abscissas and
weights are determined by the gauleg routine given by Press90.
The iterative process can now commence. Any chosen optimizing method is the main part of this process, with underlying algorithms such as integration, linear searches, backtracking, and matrix inversion--depending on the choice of optimizing method. Two linear searches were tested: Davidon's method Walsh, which is a cubic interpolation method, and the secant method Shew, which is a version of the Newton-Raphson method in one dimension where the second derivative has been eliminated. There was practically no difference in the performance between these two methods, so the secant method was applied, mainly due to its simplicity. The source code for the secant method can be found in Appendix A.1
The linear searches include a number of calculations of both the action and its gradient
at each iteration, which in some cases can take up a considerable amount of computing
time. Therefore, the linear search is sometimes omitted, simply using a preassigned
value for the step length factor 
(see equation 4.5), which either
is gradually set to unity when closing in on the solution when applied for the
Newton-like methods, or is just modified so that the action decreases at each iteration
(Steepest Descent method). This may result in more iterations, but the total computing
time can be reduced. The backtracking scheme is, in reality, just a more advanced version
of these rather crude stepsize algorithms, where the step length is determined to be
equal to or less than unity by fulfilling certain criteria for the decrease of the
action (see equation 5.21). The computational effort of this method is
less than for the linear search, but somewhat more than for the crude stepsize
algorithms, and it only applies for the Newton-Raphson and the Secant method.
The Newton-Raphson optimizing method needs a matrix inversion algorithm; i.e. solving a linear system as shown in equation (4.18). There are several methods for this type of problem, e.g. Gauss-Jordan elimination, Gauss elimination with backsubstitution, LU decomposition, and Cholesky factorization. All these four methods have been tried here, and not surprisingly, since the matrix inverse are not needed per se, the LU decomposition and Cholesky factorization proved to be the most efficient. The Cholesky factorization is the better of the two, but has the limitation that the matrix must be positive definite. Therefore, the Cholesky factorization will only be used when the Hessian is positive definite, otherwise, the LU decomposition will be applied. More on these methods in Press.
The iterative process was terminated when the minimum point (or saddle point) is reached
within a given accuracy or, optionally, when a maximum number of iterations is exceeded.
A measure for the approach to a stationary point is when the norm of the gradient becomes
sufficiently small, that is, when a point 
is obtained with
where 
is a small positive scalar. If the solution is to be
determined to machine accuracy in double precision, the iterations are terminated when
, which means that the average error for each variable is
about 
. This accuracy is somewhat extravagant, since the error in the input
parameters and the error due to the truncation in the parameterization of the orbits in
equation (3.21) is at least ten factor of ten larger. But, for some of the
optimizing methods, e.g. the Newton-Raphson method, there is negligible extra
computational cost in driving the solution to machine accuracy; when the accuracy of
about 
has been reached, the method converge to machine accuracy in just a few
iterations. For other optimizing methods with slower convergence, e.g. Steepest Descent,
a more modest accuracy must be accepted in order to avoid impossibly long computing
time.
After the termination of the iterative process, the code continues to calculate the orbits and physical velocities. In plotting the orbits, we are only interested in the comoving positions, which are given by (3.21); i.e.
The trial functions 
can now be calculated at arbitrary times in the
expansion parameter since there now is no need for any numerical integration. 500
equally spaced points in time were used here, which results in a very smooth curve.
As a test of the validity of the AVP orbits, a conventional N-body code is often used.
Contrary to the AVP, these codes solve initial value problems, where both the velocity
and position of each particle at a given point in time is needed. The N-body codes then
calculate the orbits forward in time. Therefore, to be able to check the AVP orbits
against a N-body program, we need both the position and velocity at a given time early
in each mass tracer's history from the former as input to the latter, then in physical
coordinates.
(More on N-body codes
in Section 5.2.) The physical positions are given by equation
(2.5) and the physical velocities by equation (2.6).The
in the latter is either given by equation (3.33) or
(3.34). The peculiar velocity given in the last term on the right hand
side of equation (2.6) requires that the time derivatives of the comoving
positions are determined. By equation (5.5), we have
which enable us to determine the physical velocity at any given time. The
physical velocities at the current time are of special interest, since they can be
compared to observed radial velocities. For example, the observed value of the radial
velocity of M31 can be used to determine its mass, and thereby determine the mass of all
the other mass tracers in the system, since they are scaled to the mass of M31. This is
one of the two procedures used to determine the mass of the system: the mass of M31 is
adjusted until there is an acceptable fit between the observed and the predicted radial
velocity of M31, which is -123 km s 
P-89. In the second procedure, one
tries to get an best overall fit between observed and predicted radial
velocities for several mass tracers by adjusting the mass of M31, using mass tracers
that we have acceptably accurate observations of.
In this section, the orbits are intended to be compared with the solutions of P-89, which means that the coordinate system has to be rotated, so that the connecting line between the MW and M31 is aligned with the x-axis. Hence, both the positions and the velocities are adjusted to this transformed coordinate system. More on the presentation of the orbits in the next subsection.
After this slight adjustment, the output routine is being reached. Firstly, the number of
iterations and the norm of the gradients are written to screen. Secondly, the
type of solution arrived at is being determined: minimum, maximum. or saddle point. This
is done by looking at the eigenvalues of the Hessian matrix at the stationary point: if
the eigenvalues are all positive, we have a minimum; if they are all positive, we have a
maximum; and if there are both positive and negative eigenvalues, we are dealing with a
saddle point. The eigenvalues are calculated by using the tred2 and tqli from
Press90, where the former reduce the symmetric matrix to a tridiagonal system,
while the latter determine the eigenvalues of such a tridiagonal matrix.
The closing stages of the AVP code has now been reached, were the results are written to file: the coefficients, the orbits (that is, the positions at 500 points in time), and the initial values used as input to the N-body code.
I close this section by recapitulating the expressions for the action and its derivative
with respect to the coefficients :
and
where
Notice that this is the derivative prior to the integration by parts in
equation (3.24). We also need an expression for the second
derivative, used by the Newton-Raphson optimizing method. Derivation of
(5.8) with respect to all the other coefficients 
gives:
where
with 
, and being the common
Kronecker delta. Note the difference between the subscript m and the mass 
.
In the equations above, the cutoff length c discussed in Section 3.1 is not included, but is being used in the calculations of the orbits. The purpose of the cutoff is to avoid singularities in the gravitational acceleration at zero separation, and to simulate spherical halos. The cutoff, often also referred to as a softening parameter, is introduced to the first term in equation (5.9) in the following way:
The second term in equation (5.7) is being modified in the
same manner. In this thesis the cutoff of P-90 is adapted; i.e. 
of the
radius given in equation (3.36). This is a relatively large cutoff, and
should allow reasonably large halos.
