Last change
on this file since 7837 was
7837,
checked in by mungkasi, 14 years ago

Again, adding some codes for 1d problems on debris avalanche and periodic waves.

File size:
1.1 KB

Rev  Line  

[7837]  1  ## module gaussPivot 

 2  ''' x = gaussPivot(a,b,tol=1.0e9). 

 3  Solves [a]{x} = {b} by Gauss elimination with 

 4  scaled row pivoting 

 5  ''' 

 6  from Numeric import * 

 7  import swap 

 8  import error 

 9  

 10  def gaussPivot(a,b,tol=1.0e12): 

 11  n = len(b) 

 12  

 13  # Set up scale factors 

 14  s = zeros((n),Float) 

 15  for i in range(n): 

 16  s[i] = max(abs(a[i,:])) 

 17  

 18  for k in range(0,n1): 

 19  

 20  # Row interchange, if needed 

 21  p = int(argmax(abs(a[k:n,k])/s[k:n])) + k 

 22  if abs(a[p,k]) < tol: error.err('Matrix is singular') 

 23  if p != k: 

 24  swap.swapRows(b,k,p) 

 25  swap.swapRows(s,k,p) 

 26  swap.swapRows(a,k,p) 

 27  

 28  # Elimination 

 29  for i in range(k+1,n): 

 30  if a[i,k] != 0.0: 

 31  lam = a[i,k]/a[k,k] 

 32  a[i,k+1:n] = a [i,k+1:n]  lam*a[k,k+1:n] 

 33  b[i] = b[i]  lam*b[k] 

 34  if abs(a[n1,n1]) < tol: error.err('Matrix is singular') 

 35  

 36  # Back substitution 

 37  for k in range(n1,1,1): 

 38  b[k] = (b[k]  dot(a[k,k+1:n],b[k+1:n]))/a[k,k] 

 39  return b 

 40  

 41  

 42  

 43  

Note: See
TracBrowser
for help on using the repository browser.