# "Numerical Recipes" ludcmp() C code on page 46 translated into Python # I make it as simple as possible, disregard efficiency. # Here a is a list of list, n is integer (size of the matrix) # index is a list, and d is also a list of size 1 # Python list index starts from 0. So matrix index is from 0 to n-1. # import math def ludcmp(a, n, indx, d): d[0] = 1.0 # looking for the largest a in each row and store it in vv as inverse # We need a new list same size as indx, for this we use .copy() vv=indx.copy() for i in range(0,n): big = 0.0 for j in range(0,n): temp = math.fabs(a[i][j]) if (temp > big): big = temp vv[i]=1.0/big # # run Crout's algorithm for j in range(0,n): # top half & bottom part are combined # but the upper limit l for k sum is different big = 0.0 for i in range(0,n): if(i=j): dum = vv[i]*math.fabs(sum) if (dum >= big): big = dum imax = i # end if (i>= ...) # end for i # pivoting part, swap row j with row imax, a[j] is a whole row if (j != imax): dum = a[imax] a[imax] = a[j] a[j] = dum d[0] = - d[0] vv[imax] = vv[j] # end if (j != ...) # divide by the beta diagonal value indx[j] = imax dum = 1.0/a[j][j] for i in range(j+1,n): a[i][j] *= dum # end for i # end for j # end of def ludcmp # We do backward substitution in lubksb() take the row swapped LU decomposed # a, size n, and swapping indx, and b vector as input. The output is # in b after calling. def lubksb(a, n, indx, b): ii = -1 # forward for i in range(0,n): ip = indx[i] sum = b[ip] b[ip] = b[i] if(ii != -1): for j in range(ii,i): sum -= a[i][j]*b[j] elif(sum != 0): ii = i b[i] = sum # bote alpha_{ii} is 1 above # backward for i in range(n-1,-1,-1): sum = b[i] for j in range(i+1,n): sum -= a[i][j]*b[j] b[i] = sum/a[i][i] # end lubksb() # unfortunately a is destroyed (become swapped LU) def linearsolver(a,n,b): indx = list(range(n)) d =[1] ludcmp(a,n,indx,d) x = b.copy() lubksb(a,n,indx,x) print("x=",x) return x # end linearsolver a=[[1,3,3,-5],[2,-4,7,-1],[7,0.5,3,-6],[9,-2,3,8]] n=4 b=[0,2,3,-10] linearsolver(a,n,b)