So I have a problem with inverting matrices.
I'm trying to invert a 3x3 matrix via gaussian elimination, and use this to solve a system of linear equations. I have a python script to do so, that's called as part of a class. I've tried using arbitrary precision via the mpmath module, but the problem persists. The matrix that I want to invery is:
14.0,2746627.0,170.6155733706259951532047125510871410369873046875
2746627.0,538854462769.0,33472668.1214660549849622128704140777699649333953857421875
170.6155733706259951532047125510871410369873046875,33472668.1214660549849622128704140777699649333953857421875,2079.2624245838425583639446197230982400395888676326095937481
And my code is:
def invertmat(self,locstring):
"""Here we will invert the matrix that we have found using RRE"""
if usempf==1:
self.invmat=[[mp.mpf(0) for i in range(veclen)] for j in range(veclen)]
for i in range(veclen):
self.invmat[i][i]=mp.mpf(1)
else:
self.invmat=[[0.0 for i in range(veclen)] for j in range(veclen)]
for i in range(veclen):
self.invmat[i][i]=1
if declarestages==1: print "Inverting matrix..."
#we will use RRE, so we will put the first column in RRE form, then the second, etc
#this means that we will iterate over the column and then the row
for j in range(veclen):
for i in range(veclen):
#first we have to check whether or not the term on the leading diagonal is zero
#if it is, we zero out the row and column
if self.mat[j][j]==0:
if usempf==1:
self.mat[i][j]=mp.mpf(0.0)
self.invmat[i][j]=mp.mpf(0.0)
self.vec[j]=mp.mpf(0.0)
self.badvec[j]=mp.mpf(0.0)
else:
self.mat[i][j]=float(0.0)
self.invmat[i][j]=float(0.0)
self.vec[j]=float(0.0)
self.badvec[j]=float(0.0)
else:
if i==j:
#in this case, we are on the leading diagonal
#so we will make the leading diagonal zero
self.vec[i]=self.vec[i]/self.mat[i][i]
self.badvec[i]=self.badvec[i]/self.mat[i][i]
print "Dividing by "+str(self.mat[i][i])
for k in range(veclen):
if k!=j:
self.invmat[i][k]=self.invmat[i][k]/self.mat[i][i]
self.mat[i][k]=self.mat[i][k]/self.mat[i][i]
self.invmat[i][i]=self.invmat[i][i]/self.mat[i][i]
self.mat[i][i]=self.mat[i][i]/self.mat[i][i]
if i!=j:
#in this case, we are not on the leading diagonal
#so we want this element to be zero
#we will use the fact that the leading diagonal is nonzero
#we only care about rounding to zero on the leading diagonal, as that's the only thing we divide by
mult=self.mat[i][j]/self.mat[j][j]
self.vec[i]-=self.vec[j]*mult
self.badvec[i]-=self.badvec[j]*mult
for k in range(veclen):
if k!=j:
self.mat[i][k]-=self.mat[j][k]*mult
self.invmat[i][k]-=self.invmat[j][k]*mult
self.mat[i][j]-=self.mat[j][j]*mult
self.invmat[i][j]-=self.invmat[j][j]*mult
if writeprogresstofile>0:
f.write("vec"+str(j)+":\n")
for i in self.vec:
f.write(str(i)+",")
f.write("\nmat"+str(j)+":\n")
for i in self.mat:
for k in i:
f.write(str(k)+",")
f.write("\n")
f.write("invmat"+str(j)+":\n")
for i in self.invmat:
for k in i:
f.write(str(k)+",")
f.write("\n")
the relevant variables are that veclen is 3, and f is a file that I write progress to.
When I invert this, however, I get:
271635475736322.3449137725807845490330514652153613387246069,123760956.39505916619461152377801645592270193200911064532913,-24281616049723.701885662857417315428214747782727699841050359
123760956.39505916619461152377801645592270193200911064830143,56.387244784681245388662458884703294053440414360396328263187,-11063047.066704691682385396698842928335615664887677997376554
-24281616049723.701885662857417315428214747782727699841098208,-11063047.06670469168238539669884292833561566488767799713266,2170544465289.5403113558612668273593392823016506598570495833
which is not correct, as multiplying this by the original gives:
-4,1.3113E-05,0.25
393216,4.25,-81920
-8,0.000190735,-1
What is the source of my problem?