I'm attempting to simulate a certain physical system. In order to propagate solutions I need to be able to multiply matrices of determinant = 1 which describe each part of the system. In the code below T(variables) is a 2-dimensional matrix with det(T) = 1, i just indicates the region number and the rest is irrelevant.
When I run this code for systems with more than 30 regions, the final Msys no longer has determinant = 1. I checked the value of the determinant of Msys throughout the calculation and it's 1 for the first few iterations but then it starts diverging from that. I've tried putting dtype = float64 when creating the array T to see if this would improve the precision and stop it from breaking down but I saw no improvement.
Is there any way I could write the code to avoid the error accumulating or a way I can increase the amount of decimal places numpy stores so to make the error negligible for systems with 100+ regions.
for i in range(n):
if i == 0:
Msys = T(L[i],i,k)
else:
Msys = numpy.dot(T(L[i]-L[i-1],i,k), Msys)
return Msys