Thomas algorithm

This is a function that solves a triangular matrix using Thomas Algorithm. This function will be very useful in the future. I'm planing to use it to solve partial differential equations later on.

The Code:

Save to file in a file named: ChemE.py

def thomas(a,b,c,d):
 '''Uses Thomas algorithm for solving a tridiagonal matrix for n unknowns.
  a, b, and  are a list of the matrix entries
  Matrix form of:
  [b1 c1    ] [x1] [d1]
  |a2 b2 c2   | |x2|  |d2|
  | a3 b3 c3  | |x3|  |d3|
  |     | |' |= | '|
  |     | |' |  | '|
  [   an bn cn] |xn] [dn]
 '''
 # Finding the size of the matrix and determine n
 n = len(b)
 #print n # Used for debugging
 # Test the size of a and c
 if len(a) != n-1:
  print 'Wrong index size for a.\n A should have an index of', n-1, '\n Your a has ', len(a)
  exit()
 if len(c) != n-1:
  print 'Wrong index size for c.\n C should have an index of', n-1, '\n Your c has', len(c)
  exit()
 
 # Converting to float and appending 0.0 to c
 for i in range(0,len(a)):
  a[i] = float(a[i])
 for i in range(0,len(b)):
  b[i] = float(b[i])
 for i in range(0,len(c)):
  c[i] = float(c[i])
 for i in range(0,len(d)):
  d[i] = float(d[i])
 c.append(0.0) # Hack to make the function to work

 # Calculate p and q
 p = []; q= [] 
 p.append(c[0]/b[0]); q.append(d[0]/b[0])
 for j in range(1,n):
  pj = c[j]/(b[j] - a[j-1]* p[j-1])
  qj = (d[j] - a[j-1]*q[j-1])/(b[j] - a[j-1]* p[j-1])
  p.append(pj); q.append(qj)
 #print p,q # Used for debugging the code!

 # Back sub
 x = []; x.append(q[n-1])
 for j in range(n-2,-1,-1):
  xj = q[j] - p[j]*x[0]  # Value holder
  x.insert(0,xj)   # Building the list backwards

 # Return the value
 return x

Example code execution:

This example is take from SEPARATION PROCESS PRINCIPLES: Chemical and Biochemical Operations 3ed by Seader, Henley, and Roper. Problem number is 10.5 on page 406.

"Use the Thomas algorithm to solve the following tridiagonal-matrix equation for the x vector."
| -6 3 0  0 0  0 0 0| |x1|    |   0   |
|  3 -4.5 3 0  0 0 0| |x2|    |   0   |
|  0  1.5 -7.5 3 0 0| |x3| = | 100 |
|  0  0  4.5 -7.5   3| |x4|    |   0   |
|  0   0  0   4.5-4.5| |x5|    |   0   |

#!/usr/bin/python

# Import the function:
from ChemE import thomas
# Matrix
a= [3,1.5,4.5,4.5]
b= [-6,-4.5,-7.5,-7.5,-4.5] 
c= [3,3,3,3]
d= [0,0,100,0,0]
# Execute the function
x = thomas(a,b,c,d)
# Print the result
print x

First post and first code. Primefinder

The first post on this blog. I have previously used math packages, namely MATLAB, to solve chemical engineering problems. In an attempt to move away from commercial products and gain a better understanding of what is going on behind the curtain; I've decided to learn how to use python to solve these problems. I'm posting the code as I'm learning. In the beginning it will be mostly python examples but I hope I soon can get in to solve real problems.

This code example generates all prime numbers within a given range. It is a good python example with uses: if statements, for loops, while loops, and user inputs. I know the code can be optimized but I leave it to you. The code is licensed under GNU GPLv3.

Copy the code and save it to a file named primefinder.py and run.

The Code:

#!/usr/bin/env python

print 'This program finds all prime numbers between a range given'
nMin = eval(raw_input('Enter left hand side of the range: '))
nMax = eval(raw_input('Enter right hand side, number is excluded: '))

# Checking if nMin is to small for this algorithm
if nMin <= 2: # If less than 2, add to and start with 3
 a = [2] 
 print '2' 
 nMin = 3
else:
 a = []

i = 2
# Loop to find all primes in the range
for n in range(nMin,nMax):
 while i < n:   # For every number less than n:
  if n%i != 0:  # If n is not a factor of i
   prime = 1 # Set the checking variable to true
   i = i +1  # Go to next number
  else: 
   prime = 0 # Otherwise set checking variable to false
   i = n + 1 # Exit the loop
 i = 2     # Reset i for next number in the range
 # Print the primes found and save them in variable a 
 if prime == 1:
  print n 
  a.append(n)

# Find out how many was found and print result
Size = len(a)
print 'The number of prime numbers found in this range is', Size

# Ask user what to do next.
while 1>0:
 print 'Would you like me to reprint all the numbers?'
 ANS = raw_input('yes/no: ')
 if ANS == 'yes':
  print a
  exit()
 elif ANS == 'no':
  exit()
 else:
  print 'Type yes or no'

An example of the code running:

This program finds all prime numbers between a range given
Enter left hand side of the range: 0
Enter right hand side, number is excluded: 50
2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
The number of prime numbers found in this range is 15
Would you like me to reprint all the numbers?
yes/no: yes
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]