solving the closed orbit fix point in pyat and AT matlab

[oscarxblanco]

Dear all,
there are some differences between the algorithm used in pyat and the one in AT matlab to find the fix point closed orbit.
In pyat the solution is found through the numpy linalg.solve function that according to the numpy documentation uses the LAPACK routine _gesv

at/pyat/at/physics/orbit.py

Lines 59 to 60 in 73e6151

	b_over_a = numpy.linalg.solve(a, b)
	r_next = ref_in - numpy.append(b_over_a, numpy.zeros((2,)))

The same step in AT matlab is done by matrix operations, and no library is called

at/atmat/atphysics/Orbit/xorbit_6.m

Line 47 in 73e6151

Ri_next = Ri + (eye(6)-J6)\(Rf-Ri-theta);

This typically leads to small differences, in most cases too small to be relevant. However, those differences could accumulate and produce a visible change in heavy computations like momentum aperture, dynamic aperture, iterative orbit corrections, long term tracking and so on.

I wonder on the importance to get the same numerical result in pyat and AT matlab, which would mean to choose one of the two methods and implement it on both versions of the toolbox.

At this moment I'm not asking or proposing modifications, I wanted to share a detail on a basic function that actually produces similar but not exactly the same numerical results in pyat and AT matlab.

[lfarv]

@oscarxblanco 👍
I know the problem, and this is why in the test sequence, in test_physics.py, the check on chromaticities has to be quite large (2.0e-4 in absolute difference). Same in atmat/attests/pytests.m.

In Matlab, the matrix "\" operator, for non-square matrices, is supposed to give the least-square solution, but with a non-documented method. From Matlab help: "A\B returns a least-squares solution to the system of equations A*x= B. x may not be the minimum-norm solution."

So the Matlab solution is difficult to implement in python. And having Matlab call directly LAPACK is also not obvious, though Matlab is supposed to be based on it…