ODE - III

With the changes in the Matrix class, another round of changes to the ODE solver. The Runge Kutta function has been primarily changed to remove all object constructors and object returns in the functions. The new solver functions is:

	def mat_ode(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, runge_vectors):
	""" Solves the ODE using Runge Kutta 4th order method. """
	mat_ode_reduce(matrix_e, matrix_a, matrix_b)
	def runge_function(x_plus_dxdt, ddt_mat, dbydt_order):
	""" Defines the function dx/dt=f(x) """
	#Calculate d/dt vector in reverse.
	for c1 in range(matrix_e.rows-1, -1, -1):
	if matrix_e.data[c1][c1]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=matrix_b.data[c1][0]*input_vector
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=0.0 # Added on Dec. 29, 2012.
	for c2 in range(matrix_b.columns):
	ddt_mat[dbydt_order-1].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	if (dbydt_order==2):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0]
	if (dbydt_order==3):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0]
	if (dbydt_order==4):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*ddt_mat[dbydt_order-2].data[c3][0]
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*x_plus_dxdt.data[c3][0]
	for c3 in range(c1+1, matrix_e.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_e.data[c1][c3]*runge_vectors[4].data[c3][0]/dt
	ddt_mat[dbydt_order-1].data[c1][0]=ddt_mat[dbydt_order-1].data[c1][0]*dt
	else:
	if not (matrix_a.data[c1][c1]):
	#The variable has no dynamics and can't even be calculated statically!
	print "This result will be wrong!"
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=0.0
	state_vector[0].data[c1][0]=0.0
	state_vector[1].data[c1][0]=0.0
	for c2 in range(matrix_a.columns):
	if not (c1==c2):
	state_vector[0].data[c1][0]+=matrix_a.data[c1][c2]*state_vector[0].data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	state_vector[0].data[c1][0]+=matrix_b.data[c1][0]*input_vector
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	state_vector[0].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	state_vector[0].data[c1][0]=-state_vector[0].data[c1][0]/matrix_a.data[c1][c1]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	return
	# Runge-Kutta 4th order method
	runge_function(state_vector[0], runge_vectors, 1) #k1=dt*(dx/dt)
	runge_function(state_vector[0], runge_vectors, 2) #k2=dt*d/dt(x+k1/2)
	runge_function(state_vector[0], runge_vectors, 3) #k3=dt*d/dt(x+k2/2)
	runge_function(state_vector[0], runge_vectors, 4) #k4=dt*d/dt(x+k3)
	for c1 in range(state_vector[1].rows):
	runge_vectors[4].data[c1][0]=(runge_vectors[0].data[c1][0]+2.0*runge_vectors[1].data[c1][0]+2.0*runge_vectors[2].data[c1][0]+runge_vectors[3].data[c1][0])*(1/6.0)
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+runge_vectors[4].data[c1][0]
	return

view raw solver_func.py hosted with ❤ by GitHub

The full code is:

	#! /usr/bin/env python
	import sys, math, matrix
	def row_op(matrix_a, matrix_b):
	"""Function reduces a set of a set of equations Ax=b to upper triangular form."""
	# Check if it is a genuine set of equations
	if not (matrix_a.rows==matrix_b.rows):
	print "Matrices do not have the same number of rows."
	else:
	for c1 in range(matrix_a.rows):
	#Check if the diagonal element is zero.
	if not matrix_a.data[c1][c1]:
	#If diagonal element is zero try to interchange it with a lower row.
	for c2 in range(c1+1, matrix_a.rows):
	if matrix_a.data[c2][c1]:
	for c3 in range(matrix_a.columns):
	matrix_a.data[c1][c3], matrix_a.data[c2][c3] = matrix_a.data[c2][c3], matrix_a.data[c1][c3]
	for c3 in range(matrix_b.columns):
	matrix_b.data[c1][c3], matrix_b.data[c2][c3] = matrix_b.data[c2][c3], matrix_b.data[c1][c3]
	#Make the diagonal element = 1.
	div_row=matrix_a.data[c1][c1]
	#Even after all possible row manipulations,
	#a diagonal element is zero, it means rank<size of matrix.
	if not div_row:
	print "Infinite solutions possible for this set of equations."
	else:
	for c2 in range(matrix_a.columns):
	matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/div_row
	for c2 in range(matrix_b.columns):
	matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/div_row
	#Row operations to remove row entries below the diagonal unity element.
	for c2 in range(c1+1, matrix_a.rows):
	mul_column=matrix_a.data[c2][c1]
	for c3 in range(matrix_a.columns):
	matrix_a.data[c2][c3]=matrix_a.data[c2][c3]-matrix_a.data[c1][c3]*mul_column
	for c3 in range(matrix_b.columns):
	matrix_b.data[c2][c3]=matrix_b.data[c2][c3]-matrix_b.data[c1][c3]*mul_column
	return
	def mat_equation(matrix_a, matrix_b, inp_signal, matrix_x):
	""" Function returns result x for equations Ax=bu """
	row_op(matrix_a, matrix_b)
	#Calculate result matrix in reverse.
	for c1 in range(matrix_a.rows-1, -1, -1):
	if matrix_a.data[c1][c1]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	matrix_x.data[c1][0]=matrix_b.data[c1][0]*input_vector
	else:
	if not (matrix_b.columns==input_signal.rows):
	print "Dimension of input vector incorrect."
	else:
	matrix_x.data[c1][0]=0.0
	for c2 in range(matrix_b.columns):
	matrix_x.data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	for c3 in range(c1+1, matrix_a.columns):
	matrix_x.data[c1][0]-=matrix_a.data[c1][c3]*matrix_x.data[c3][0]
	else:
	#If rank<size, set the (size-rank) solutions to zero.
	matrix_x.data[c1][0]=0.0
	return
	def mat_ode_reduce(matrix_e, matrix_a, matrix_b):
	""" Function to modify the matrices for E d/dt(x)=Ax+bu.
	Upper triangularization is used to calculate the value of d/dt(x)."""
	# Check if it is a genuine set of equations
	if not ((matrix_e.rows==matrix_a.rows) and (matrix_a.rows==matrix_b.rows)):
	print "Matrices do not have the same number of rows."
	else:
	for c1 in range(matrix_e.rows):
	#Check if the diagonal element is zero.
	if not matrix_e.data[c1][c1]:
	#If diagonal element is zero try to interchange it with a lower row.
	for c2 in range(c1+1, matrix_e.rows):
	if matrix_e.data[c2][c1]:
	for c3 in range(matrix_e.columns):
	matrix_e.data[c1][c3], matrix_e.data[c2][c3] = matrix_e.data[c2][c3], matrix_e.data[c1][c3]
	for c3 in range(matrix_a.columns):
	matrix_a.data[c1][c3], matrix_a.data[c2][c3] = matrix_a.data[c2][c3], matrix_a.data[c1][c3]
	for c3 in range(matrix_b.columns):
	matrix_b.data[c1][c3], matrix_b.data[c2][c3] = matrix_b.data[c2][c3], matrix_b.data[c1][c3]
	#Make the diagonal element = 1.
	div_row=matrix_e.data[c1][c1]
	#Even after all possible row manipulations,
	#a diagonal element is zero, it means rank<size of matrix.
	if not div_row:
	# print "Check the circuit. Some variables have no dynamics."
	pass
	else:
	for c2 in range(matrix_e.columns):
	matrix_e.data[c1][c2]=matrix_e.data[c1][c2]/div_row
	for c2 in range(matrix_a.columns):
	matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/div_row
	for c2 in range(matrix_b.columns):
	matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/div_row
	#Row operations to remove row entries below the diagonal unity element.
	for c2 in range(c1+1, matrix_e.rows):
	mul_column=matrix_e.data[c2][c1]
	for c3 in range(matrix_e.columns):
	matrix_e.data[c2][c3]=matrix_e.data[c2][c3]-matrix_e.data[c1][c3]*mul_column
	for c3 in range(matrix_a.columns):
	matrix_a.data[c2][c3]=matrix_a.data[c2][c3]-matrix_a.data[c1][c3]*mul_column
	for c3 in range(matrix_b.columns):
	matrix_b.data[c2][c3]=matrix_b.data[c2][c3]-matrix_b.data[c1][c3]*mul_column
	return
	def mat_ode(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, runge_vectors):
	""" Solves the ODE using Runge Kutta 4th order method. """
	mat_ode_reduce(matrix_e, matrix_a, matrix_b)
	def runge_function(x_plus_dxdt, ddt_mat, dbydt_order):
	""" Defines the function dx/dt=f(x) """
	#Calculate d/dt vector in reverse.
	for c1 in range(matrix_e.rows-1, -1, -1):
	if matrix_e.data[c1][c1]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=matrix_b.data[c1][0]*input_vector
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=0.0 # Added on Dec. 29, 2012.
	for c2 in range(matrix_b.columns):
	ddt_mat[dbydt_order-1].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	if (dbydt_order==2):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0]
	if (dbydt_order==3):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0]
	if (dbydt_order==4):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*ddt_mat[dbydt_order-2].data[c3][0]
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*x_plus_dxdt.data[c3][0]
	for c3 in range(c1+1, matrix_e.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_e.data[c1][c3]*runge_vectors[4].data[c3][0]/dt
	ddt_mat[dbydt_order-1].data[c1][0]=ddt_mat[dbydt_order-1].data[c1][0]*dt
	else:
	if not (matrix_a.data[c1][c1]):
	#The variable has no dynamics and can't even be calculated statically!
	print "This result will be wrong!"
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=0.0
	state_vector[0].data[c1][0]=0.0
	state_vector[1].data[c1][0]=0.0
	for c2 in range(matrix_a.columns):
	if not (c1==c2):
	state_vector[0].data[c1][0]+=matrix_a.data[c1][c2]*state_vector[0].data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	state_vector[0].data[c1][0]+=matrix_b.data[c1][0]*input_vector
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	state_vector[0].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	state_vector[0].data[c1][0]=-state_vector[0].data[c1][0]/matrix_a.data[c1][c1]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	return
	# Runge-Kutta 4th order method
	runge_function(state_vector[0], runge_vectors, 1) #k1=dt*(dx/dt)
	runge_function(state_vector[0], runge_vectors, 2) #k2=dt*d/dt(x+k1/2)
	runge_function(state_vector[0], runge_vectors, 3) #k3=dt*d/dt(x+k2/2)
	runge_function(state_vector[0], runge_vectors, 4) #k4=dt*d/dt(x+k3)
	for c1 in range(state_vector[1].rows):
	runge_vectors[4].data[c1][0]=(runge_vectors[0].data[c1][0]+2.0*runge_vectors[1].data[c1][0]+2.0*runge_vectors[2].data[c1][0]+runge_vectors[3].data[c1][0])*(1/6.0)
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+runge_vectors[4].data[c1][0]
	return
	R1=10.0
	R2=10.0
	R3=10.0
	L1=20.0e-3
	L2=0.0e-3
	L3=0.0e-3
	a=matrix.Matrix(2,2)
	b=matrix.Matrix(2)
	b.data=[[1.0],[0.0]]
	a.data=[[R1+R2, -R2], [-R2, R2+R3]]
	e=matrix.Matrix(2,2)
	e.data=[[L1+L2, -L2], [-L2, L2+L3]]
	curr_state_vec = matrix.Matrix(2)
	curr_state_vec.zeros(2)
	next_state_vec = matrix.Matrix(2)
	next_state_vec.zeros(2)
	dt=20.0e-6
	f=open("data_file.dat","w")
	input_volt=matrix.Matrix(2)
	ode_k1=matrix.Matrix(2)
	ode_k2=matrix.Matrix(2)
	ode_k3=matrix.Matrix(2)
	ode_k4=matrix.Matrix(2)
	ode_dbydt=matrix.Matrix(2)
	ode_solver=[ode_k1, ode_k2, ode_k3, ode_k4, ode_dbydt]
	for time_step in range(0,10000):
	t=time_step*dt
	u=100.0*math.sin(2*60*math.pi*t)
	mat_ode(e, a, b, [curr_state_vec, next_state_vec], u, dt, ode_solver)
	curr_state_vec=next_state_vec
	f.write("%s \t %s \t %s \t %s \n" %(str(t), str(u), str(next_state_vec.data[0][0]), str(next_state_vec.data[1][0])))
	print a
	print e
	print b
	f.close()

view raw solver.py hosted with ❤ by GitHub

Matrices - II

A couple of changes to the matrix.py module. The ODE solver gave me the idea of how important it was to cut down on copy constructors and destructors. So there were two functions in the Matrix class that would come under that category and will be used almost everywhere.

	#! /usr/bin/env python
	import sys
	class Matrix:
	""" For defining arrays of any given size.
	Contains overloaded operators for print, add, subtract and multiply. Does not contain = assignment operator! Use a(b) instead of a=b.
	Last updated - Dec 30, 2012 """
	def resize(self, rows, columns=1):
	""" Changes the dimension of an existing matrix by add/deleting difference of rows and columns."""
	if (rows<self.rows):
	for c1 in range(self.rows-1,rows-1,-1):
	del self.data[c1]
	self.rows=rows
	if (columns<self.columns):
	for c1 in range(self.rows):
	for c2 in range(self.columns-1, columns-1, -1):
	del self.data[c1][c2]
	for c1 in range(self.rows, rows):
	row_vector=[]
	for c2 in range(columns):
	row_vector.append(0)
	self.data.append(row_vector)
	else:
	for c1 in range(self.rows):
	for c2 in range(self.columns, columns):
	self.data[c1].append(0)
	for c1 in range(self.rows, rows):
	row_vector=[]
	for c2 in range(columns):
	row_vector.append(0)
	self.data.append(row_vector)
	self.rows=rows
	self.columns=columns
	return
	def zeros(self, rows, columns=1):
	""" Creates a zero matrix
	Dimensions have to be passed as <Matrix>.zeros(rows, columns).
	This function also creates the actual array. So every method of the class/object should call this first before using the array."""
	try:
	if self.data: pass
	except:
	self.data=[]
	for count1 in range(self.rows):
	row_vector=[]
	for count2 in range(self.columns):
	row_vector.append(0.0)
	self.data.append(row_vector)
	if not (self.rows==rows and self.columns==columns):
	self.resize(rows, columns)
	# self.data=[] Removed on Dec. 30 due to resize() function addition to class.
	for count1 in range(self.rows):
	#row_vector=[] Removed on Dec. 30 due to resize() function addition to class.
	for count2 in range(self.columns):
	#row_vector.append(0.0) Removed on Dec. 30 due to resize() function addition to class.
	self.data[count1][count2]=0.0
	#self.data.append(row_vector) Removed on Dec. 30 due to resize() function addition to class.
	return
	def __init__(self, rows=1, columns=1):
	"""Creates an object with a zero array of size rowsXcolumns.
	If only one dimension is specified specified, a column vector is assumed. If no dimensions are specified, an null matrix is created."""
	self.rows=1
	self.columns=1
	self.zeros(rows, columns)
	def __call__(self, matrix1):
	""" Replaces the assignment operator =. Use as c(b) instead of c=b.
	The = will only result in a reference to the object and not an actual copy. This method ensures it. """
	#self.rows=matrix1.rows Removed on Dec. 30 due to resize() function addition to class.
	#self.columns=matrix1.columns Removed on Dec. 30 due to resize() function addition to class.
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	self.resize(matrix1.rows, matrix1.columns)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	self.data[count1][count2]=matrix1.data[count1][count2]
	def identity(self, rows):
	""" Generates an identity matrix of rowsXrows. """
	self.rows=rows
	self.columns=rows
	self.zeros(self.rows, self.columns)
	for count in range(self.rows):
	self.data[count][count]=1.0
	def display(self):
	""" Explicit method for displaying an array. """
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	print self.data[count1][count2],
	print "\t",
	print "\n",
	def __repr__(self):
	""" Displays an array when print <Matrix> is used. """
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	print self.data[count1][count2],
	print "\t",
	print "\n",
	return " "
	def __add__(self, matrix1):
	""" Left addition of two matrices. """
	add_result=Matrix(self.rows, self.columns)
	try:
	if matrix1.rows: pass
	except AttributeError:
	print "*"*40
	print "Can't add a float or integer to a matrix directly."
	print "*"*40
	sys.exit("Logic Error")
	else:
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	print "Matrices not compatible for summation."
	sys.exit("Logic Error")
	else:
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	add_result.data[count1][count2]=self.data[count1][count2]+matrix1.data[count1][count2]
	return add_result
	def __radd__(self, matrix1):
	""" Right addition of two matrices. """
	add_result=Matrix(self.rows, self.columns)
	try:
	if matrix1.rows: pass
	except AttributeError:
	print "*"*40
	print "Can't add a float or integer to a matrix directly."
	print "*"*40
	sys.exit("Logic Error")
	else:
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	print "Matrices not compatible for summation."
	else:
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	add_result.data[count1][count2]=self.data[count1][count2]+matrix1.data[count1][count2]
	return add_result
	def __sub__(self, matrix1):
	""" Left subtraction of two matrices. """
	sub_result=Matrix(self.rows, self.columns)
	try:
	if matrix1.rows: pass
	except AttributeError:
	print "*"*40
	print "Can't subtract a float or integer to a matrix directly."
	print "*"*40
	sys.exit("Logic Error")
	else:
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	print "Matrices not compatible for subtraction."
	sys.exit("Logic Error")
	else:
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	sub_result.data[count1][count2]=self.data[count1][count2]-matrix1.data[count1][count2]
	return sub_result
	def __rsub__(self, matrix1):
	""" Right subtraction of two mnatrices. """
	sub_result=Matrix(self.rows, self.columns)
	try:
	if matrix1.rows: pass
	except AttributeError:
	print "*"*40
	print "Can't subtract a float or integer to a matrix directly."
	print "*"*40
	sys.exit("Logic Error")
	else:
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	print "Matrices not compatible for subtraction."
	sys.exit("Logic Error")
	else:
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	sub_result.data[count1][count2]=matrix1.data[count1][count2]-self.data[count1][count2]
	return sub_result
	def __mul__(self, matrix1):
	""" Left multiplication of two matrices.
	The multiplication is NOT done elementwise but in the conventional mathematical sense. """
	try:
	if matrix1.rows: pass
	except AttributeError:
	mul_result=Matrix(self.rows, self.columns)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	mul_result.data[count1][count2]=self.data[count1][count2]*matrix1
	else:
	mul_result=Matrix()
	if not (self.columns==matrix1.rows):
	print "Matrices not compatible for multiplication."
	else:
	mul_result.zeros(self.rows, matrix1.columns)
	for count1 in range(self.rows):
	for count2 in range(matrix1.columns):
	for count3 in range(self.columns):
	mul_result.data[count1][count2]+=self.data[count1][count3]*matrix1.data[count3][count2]
	return mul_result
	def __rmul__(self, matrix1):
	""" Right multiplication of two matrices.
	The multiplication is NOT done elementwise but in the conventional mathematical sense. """
	try:
	if matrix1.rows: pass
	except AttributeError:
	mul_result=Matrix(self.rows, self.columns)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	mul_result.data[count1][count2]=self.data[count1][count2]*matrix1
	else:
	mul_result=Matrix()
	if not (self.rows==matrix1.columns):
	print "Matrices not compatible for multiplication."
	else:
	mul_result.zeros(matrix1.rows, self.columns)
	for count1 in range(matrix1.rows):
	for count2 in range(self.columns):
	for count3 in range(self.rows):
	mul_result.data[count1][count2]+=self.data[count1][count3]*matrix1.data[count3][count2]
	return mul_result
	def transpose(self):
	""" Returns a matrix that is transpose of the calling object array. A new object is returned. Original object is untouched. """
	trans_result=Matrix(self.columns, self.rows)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	trans_result.data[count2][count1]=self.data[count1][count2]
	return trans_result

view raw matrix.py hosted with ❤ by GitHub

The function added has been the "resize" function. So if a matrix exists but is now equated to another matrix of another dimension, the object will be modified rather than deleted and built again.

	def resize(self, rows, columns=1):
	""" Changes the dimension of an existing matrix by add/deleting difference of rows and columns."""
	if (rows<self.rows):
	for c1 in range(self.rows-1,rows-1,-1):
	del self.data[c1]
	self.rows=rows
	if (columns<self.columns):
	for c1 in range(self.rows):
	for c2 in range(self.columns-1, columns-1, -1):
	del self.data[c1][c2]
	for c1 in range(self.rows, rows):
	row_vector=[]
	for c2 in range(columns):
	row_vector.append(0)
	self.data.append(row_vector)
	else:
	for c1 in range(self.rows):
	for c2 in range(self.columns, columns):
	self.data[c1].append(0)
	for c1 in range(self.rows, rows):
	row_vector=[]
	for c2 in range(columns):
	row_vector.append(0)
	self.data.append(row_vector)
	self.rows=rows
	self.columns=columns
	return

view raw mat1.py hosted with ❤ by GitHub

And the zeros function which would initialize self.data in the beginning and then rebuild it now will check whether the object exists. If it doesn't exist, it will initialize and fill it. If the new dimension is different from the existing dimensions,  it will resize it.

	def zeros(self, rows, columns=1):
	""" Creates a zero matrix
	Dimensions have to be passed as <Matrix>.zeros(rows, columns).
	This function also creates the actual array. So every method of the class/object should call this first before using the array."""
	try:
	if self.data: pass
	except:
	self.data=[]
	for count1 in range(self.rows):
	row_vector=[]
	for count2 in range(self.columns):
	row_vector.append(0.0)
	self.data.append(row_vector)
	if not (self.rows==rows and self.columns==columns):
	self.resize(rows, columns)
	# self.data=[] Removed on Dec. 30 due to resize() function addition to class.
	for count1 in range(self.rows):
	#row_vector=[] Removed on Dec. 30 due to resize() function addition to class.
	for count2 in range(self.columns):
	#row_vector.append(0.0) Removed on Dec. 30 due to resize() function addition to class.
	self.data[count1][count2]=0.0
	#self.data.append(row_vector) Removed on Dec. 30 due to resize() function addition to class.
	return

view raw mat2.py hosted with ❤ by GitHub

The __init__ function calls the function zeros with the express arguments of the actual rows and columns of the matrix different from default rows and columns so that a resize takes place.

	def __init__(self, rows=1, columns=1):
	"""Creates an object with a zero array of size rowsXcolumns.
	If only one dimension is specified specified, a column vector is assumed. If no dimensions are specified, an null matrix is created."""
	self.rows=1
	self.columns=1
	self.zeros(rows, columns)

view raw mat3.py hosted with ❤ by GitHub

Finally, the __call__ function has been changed to check if the dimensions are different and resize if so otherwise just equate them.

	def __call__(self, matrix1):
	""" Replaces the assignment operator =. Use as c(b) instead of c=b.
	The = will only result in a reference to the object and not an actual copy. This method ensures it. """
	#self.rows=matrix1.rows Removed on Dec. 30 due to resize() function addition to class.
	#self.columns=matrix1.columns Removed on Dec. 30 due to resize() function addition to class.
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	self.resize(matrix1.rows, matrix1.columns)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	self.data[count1][count2]=matrix1.data[count1][count2]

view raw mat4.py hosted with ❤ by GitHub

The addition, subtraction and multiplication operations still need a constructor for the result and the functions will then return that result. But I will try to reduce the need for those in simulations but keep them for error checking as it makes it much easier to type out these operators as regular commands.

ODE Solver - II

So my previous post wasn't the last post of 2012. The streets are full of snow and the cold is keeping me indoors with nothing else to do but code. A couple of changes to the previous solver program:

	#! /usr/bin/env python
	import sys, math, matrix
	def row_op(matrix_a, matrix_b):
	"""Function reduces a set of a set of equations Ax=b to upper triangular form."""
	# Check if it is a genuine set of equations
	if not (matrix_a.rows==matrix_b.rows):
	print "Matrices do not have the same number of rows."
	else:
	for c1 in range(matrix_a.rows):
	#Check if the diagonal element is zero.
	if not matrix_a.data[c1][c1]:
	#If diagonal element is zero try to interchange it with a lower row.
	for c2 in range(c1+1, matrix_a.rows):
	if matrix_a.data[c2][c1]:
	for c3 in range(matrix_a.columns):
	matrix_a.data[c1][c3], matrix_a.data[c2][c3] = matrix_a.data[c2][c3], matrix_a.data[c1][c3]
	for c3 in range(matrix_b.columns):
	matrix_b.data[c1][c3], matrix_b.data[c2][c3] = matrix_b.data[c2][c3], matrix_b.data[c1][c3]
	#Make the diagonal element = 1.
	div_row=matrix_a.data[c1][c1]
	#Even after all possible row manipulations,
	#a diagonal element is zero, it means rank<size of matrix.
	if not div_row:
	print "Infinite solutions possible for this set of equations."
	else:
	for c2 in range(matrix_a.columns):
	matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/div_row
	for c2 in range(matrix_b.columns):
	matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/div_row
	#Row operations to remove row entries below the diagonal unity element.
	for c2 in range(c1+1, matrix_a.rows):
	mul_column=matrix_a.data[c2][c1]
	for c3 in range(matrix_a.columns):
	matrix_a.data[c2][c3]=matrix_a.data[c2][c3]-matrix_a.data[c1][c3]*mul_column
	for c3 in range(matrix_b.columns):
	matrix_b.data[c2][c3]=matrix_b.data[c2][c3]-matrix_b.data[c1][c3]*mul_column
	return
	def mat_equation(matrix_a, matrix_b, inp_signal, matrix_x):
	""" Function returns result x for equations Ax=bu """
	row_op(matrix_a, matrix_b)
	#Calculate result matrix in reverse.
	for c1 in range(matrix_a.rows-1, -1, -1):
	if matrix_a.data[c1][c1]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	matrix_x.data[c1][0]=matrix_b.data[c1][0]*input_vector
	else:
	if not (matrix_b.columns==input_signal.rows):
	print "Dimension of input vector incorrect."
	else:
	matrix_x.data[c1][0]=0.0 # Added on Dec. 29, 2012.
	for c2 in range(matrix_b.columns):
	matrix_x.data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	for c3 in range(c1+1, matrix_a.columns):
	matrix_x.data[c1][0]-=matrix_a.data[c1][c3]*matrix_x.data[c3][0]
	else:
	#If rank<size, set the (size-rank) solutions to zero.
	matrix_x.data[c1][0]=0.0
	return
	def mat_ode_reduce(matrix_e, matrix_a, matrix_b):
	""" Function to modify the matrices for E d/dt(x)=Ax+bu.
	Upper triangularization is used to calculate the value of d/dt(x)."""
	# Check if it is a genuine set of equations
	if not ((matrix_e.rows==matrix_a.rows) and (matrix_a.rows==matrix_b.rows)):
	print "Matrices do not have the same number of rows."
	else:
	for c1 in range(matrix_e.rows):
	#Check if the diagonal element is zero.
	if not matrix_e.data[c1][c1]:
	#If diagonal element is zero try to interchange it with a lower row.
	for c2 in range(c1+1, matrix_e.rows):
	if matrix_e.data[c2][c1]:
	for c3 in range(matrix_e.columns):
	matrix_e.data[c1][c3], matrix_e.data[c2][c3] = matrix_e.data[c2][c3], matrix_e.data[c1][c3]
	for c3 in range(matrix_a.columns):
	matrix_a.data[c1][c3], matrix_a.data[c2][c3] = matrix_a.data[c2][c3], matrix_a.data[c1][c3]
	for c3 in range(matrix_b.columns):
	matrix_b.data[c1][c3], matrix_b.data[c2][c3] = matrix_b.data[c2][c3], matrix_b.data[c1][c3]
	#Make the diagonal element = 1.
	div_row=matrix_e.data[c1][c1]
	#Even after all possible row manipulations,
	#a diagonal element is zero, it means rank<size of matrix.
	if not div_row:
	# print "Check the circuit. Some variables have no dynamics."
	pass
	else:
	for c2 in range(matrix_e.columns):
	matrix_e.data[c1][c2]=matrix_e.data[c1][c2]/div_row
	for c2 in range(matrix_a.columns):
	matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/div_row
	for c2 in range(matrix_b.columns):
	matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/div_row
	#Row operations to remove row entries below the diagonal unity element.
	for c2 in range(c1+1, matrix_e.rows):
	mul_column=matrix_e.data[c2][c1]
	for c3 in range(matrix_e.columns):
	matrix_e.data[c2][c3]=matrix_e.data[c2][c3]-matrix_e.data[c1][c3]*mul_column
	for c3 in range(matrix_a.columns):
	matrix_a.data[c2][c3]=matrix_a.data[c2][c3]-matrix_a.data[c1][c3]*mul_column
	for c3 in range(matrix_b.columns):
	matrix_b.data[c2][c3]=matrix_b.data[c2][c3]-matrix_b.data[c1][c3]*mul_column
	return
	def mat_ode(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, runge_vectors):
	""" Solves the ODE using Runge Kutta 4th order method. """
	mat_ode_reduce(matrix_e, matrix_a, matrix_b)
	def runge_function(x_plus_dxdt, ddt_mat):
	""" Defines the function dx/dt=f(x) """
	#Calculate d/dt vector in reverse.
	for c1 in range(matrix_e.rows-1, -1, -1):
	if matrix_e.data[c1][c1]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	ddt_mat.data[c1][0]=matrix_b.data[c1][0]*input_vector
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	ddt_mat.data[c1][0]=0.0 # Added on Dec. 29, 2012.
	for c2 in range(matrix_b.columns):
	ddt_mat.data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	for c3 in range(matrix_a.columns):
	ddt_mat.data[c1][0]-=matrix_a.data[c1][c3]*x_plus_dxdt.data[c3][0]
	for c3 in range(c1+1, matrix_e.columns):
	ddt_mat.data[c1][0]-=matrix_e.data[c1][c3]*ddt_mat.data[c3][0]
	else:
	if not (matrix_a.data[c1][c1]):
	#The variable has no dynamics and can't even be calculated statically!
	print "This result will be wrong!"
	else:
	ddt_mat.data[c1][0]=0.0
	state_vector[0].data[c1][0]=0.0
	state_vector[1].data[c1][0]=0.0
	for c2 in range(matrix_a.columns):
	if not (c1==c2):
	state_vector[0].data[c1][0]+=matrix_a.data[c1][c2]*state_vector[0].data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	state_vector[0].data[c1][0]+=matrix_b.data[c1][0]*input_vector
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	state_vector[0].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	state_vector[0].data[c1][0]=-state_vector[0].data[c1][0]/matrix_a.data[c1][c1]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	return ddt_mat
	# Runge-Kutta 4th order method
	k1=runge_vectors[0]
	k1(runge_function(state_vector[0], runge_vectors[4])*dt) #k1=dt*(dx/dt)
	k2=runge_vectors[1]
	k2(runge_function(state_vector[0]+0.5*k1, runge_vectors[4])*dt) #k2=dt*d/dt(x+k1/2)
	k3=runge_vectors[2]
	k3(runge_function(state_vector[0]+0.5*k2, runge_vectors[4])*dt) #k3=dt*d/dt(x+k2/2)
	k4=runge_vectors[3]
	k4(runge_function(state_vector[0]+k3, runge_vectors[4])*dt) #k4=dt*d/dt(x+k3)
	for c1 in range(state_vector[1].rows):
	state_vector[1].data[c1][0]=(state_vector[0].data[c1][0]+(k1.data[c1][0]+2.0*k2.data[c1][0]+2.0*k3.data[c1][0]+k4.data[c1][0])*(1/6.0))
	return
	R1=10.0
	R2=10.0
	R3=10.0
	L1=20.0e-3
	L2=0.0e-3
	L3=0.0e-3
	a=matrix.Matrix(2,2)
	b=matrix.Matrix(2)
	b.data=[[1.0],[0.0]]
	a.data=[[R1+R2, -R2], [-R2, R2+R3]]
	e=matrix.Matrix(2,2)
	e.data=[[L1+L2, -L2], [-L2, L2+L3]]
	curr_state_vec = matrix.Matrix(2)
	curr_state_vec.zeros(2)
	next_state_vec = matrix.Matrix(2)
	next_state_vec.zeros(2)
	dt=20.0e-6
	f=open("data_file.dat","w")
	input_volt=matrix.Matrix(2)
	ode_k1=matrix.Matrix(2)
	ode_k2=matrix.Matrix(2)
	ode_k3=matrix.Matrix(2)
	ode_k4=matrix.Matrix(2)
	ode_dbydt=matrix.Matrix(2)
	ode_solver=[ode_k1, ode_k2, ode_k3, ode_k4, ode_dbydt]
	for time_step in range(0,10000):
	t=time_step*dt
	u=100.0*math.sin(2*60*math.pi*t)
	mat_ode(e, a, b, [curr_state_vec, next_state_vec], u, 20.0e-6, ode_solver)
	curr_state_vec=next_state_vec
	f.write("%s \t %s \t %s \t %s \n" %(str(t), str(u), str(next_state_vec.data[0][0]), str(next_state_vec.data[1][0])))
	print a
	print e
	print b
	f.close()

view raw ode2.py hosted with ❤ by GitHub

The changes are a couple of bug fixes. In the function "mat_equation"
matrix_x.data[c1][0]=0.0
has been added as the result vector x has to be initialized to zero if the input is a vector. Similarly, in function "mat_ode"
ddt_mat.data[c1][0]=0.0
has been added to initialize value of dx/dt if input is a vector.


A major change has been to be able to solve the ODE when matrix E is singular. In the case, for an nXn matrix E, the rank will be k<n. This means that after the row operations, the lowest (n-k) rows will be zero rows. But if the diagonal elements of matrix A of these rows are not zero, the variable could still be calculated as that what this code does:

	if not (matrix_a.data[c1][c1]):
	#The variable has no dynamics and can't even be calculated statically!
	print "This result will be wrong!"
	else:
	ddt_mat.data[c1][0]=0.0
	state_vector[0].data[c1][0]=0.0
	state_vector[1].data[c1][0]=0.0
	for c2 in range(matrix_a.columns):
	if not (c1==c2):
	state_vector[0].data[c1][0]+=matrix_a.data[c1][c2]*state_vector[0].data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	state_vector[0].data[c1][0]+=matrix_b.data[c1][0]*input_vector
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	state_vector[0].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	state_vector[0].data[c1][0]=-state_vector[0].data[c1][0]/matrix_a.data[c1][c1]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]

view raw ode3.py hosted with ❤ by GitHub

The concept here is that through the matrix A, these variables are still linked to the dynamic variables and so even if they are not directly affected by an input signal, they should be non-zero.

Solving ODEs

One of the major components in developing the circuit simulator is solving the ordinary differential equation (ODE) “E dx/dt = Ax + bu”. Moreover, the matrices need to be determined based on the network topology, values of different components and finally in the case of power electronics, the state of a device. So to begin with, a function to solve the ODE.

	#! /usr/bin/env python
	import sys, math, matrix
	def row_op(matrix_a, matrix_b):
	"""Function reduces a set of a set of equations Ax=b to upper triangular form."""
	# Check if it is a genuine set of equations
	if not (matrix_a.rows==matrix_b.rows):
	print "Matrices do not have the same number of rows."
	else:
	for c1 in range(matrix_a.rows):
	#Check if the diagonal element is zero.
	if not matrix_a.data[c1][c1]:
	#If diagonal element is zero try to interchange it with a lower row.
	for c2 in range(c1+1, matrix_a.rows):
	if matrix_a.data[c2][c1]:
	for c3 in range(matrix_a.columns):
	matrix_a.data[c1][c3], matrix_a.data[c2][c3] = matrix_a.data[c2][c3], matrix_a.data[c1][c3]
	for c3 in range(matrix_b.columns):
	matrix_b.data[c1][c3], matrix_b.data[c2][c3] = matrix_b.data[c2][c3], matrix_b.data[c1][c3]
	#Make the diagonal element = 1.
	div_row=matrix_a.data[c1][c1]
	#Even after all possible row manipulations,
	#a diagonal element is zero, it means rank<size of matrix.
	if not div_row:
	print "Infinite solutions possible for this set of equations."
	else:
	for c2 in range(matrix_a.columns):
	matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/div_row
	for c2 in range(matrix_b.columns):
	matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/div_row
	#Row operations to remove row entries below the diagonal unity element.
	for c2 in range(c1+1, matrix_a.rows):
	mul_column=matrix_a.data[c2][c1]
	for c3 in range(matrix_a.columns):
	matrix_a.data[c2][c3]=matrix_a.data[c2][c3]-matrix_a.data[c1][c3]*mul_column
	for c3 in range(matrix_b.columns):
	matrix_b.data[c2][c3]=matrix_b.data[c2][c3]-matrix_b.data[c1][c3]*mul_column
	return
	def mat_equation(matrix_a, matrix_b, inp_signal, matrix_x):
	""" Function returns result x for equations Ax=bu """
	row_op(matrix_a, matrix_b)
	#Calculate result matrix in reverse.
	for c1 in range(matrix_a.rows-1, -1, -1):
	if matrix_a.data[c1][c1]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	matrix_x.data[c1][0]=matrix_b.data[c1][0]*input_vector
	else:
	if not (matrix_b.columns==input_signal.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	matrix_x.data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	for c3 in range(c1+1, matrix_a.columns):
	matrix_x.data[c1][0]-=matrix_a.data[c1][c3]*matrix_x.data[c3][0]
	else:
	#If rank<size, set the (size-rank) solutions to zero.
	matrix_x.data[c1][0]=0.0
	return
	def mat_ode_reduce(matrix_e, matrix_a, matrix_b):
	""" Function to modify the matrices for E d/dt(x)=Ax+bu.
	Upper triangularization is used to calculate the value of d/dt(x)."""
	# Check if it is a genuine set of equations
	if not ((matrix_e.rows==matrix_a.rows) and (matrix_a.rows==matrix_b.rows)):
	print "Matrices do not have the same number of rows."
	else:
	for c1 in range(matrix_e.rows):
	#Check if the diagonal element is zero.
	if not matrix_e.data[c1][c1]:
	#If diagonal element is zero try to interchange it with a lower row.
	for c2 in range(c1+1, matrix_e.rows):
	if matrix_e.data[c2][c1]:
	for c3 in range(matrix_e.columns):
	matrix_e.data[c1][c3], matrix_e.data[c2][c3] = matrix_e.data[c2][c3], matrix_e.data[c1][c3]
	for c3 in range(matrix_a.columns):
	matrix_a.data[c1][c3], matrix_a.data[c2][c3] = matrix_a.data[c2][c3], matrix_a.data[c1][c3]
	for c3 in range(matrix_b.columns):
	matrix_b.data[c1][c3], matrix_b.data[c2][c3] = matrix_b.data[c2][c3], matrix_b.data[c1][c3]
	#Make the diagonal element = 1.
	div_row=matrix_e.data[c1][c1]
	#Even after all possible row manipulations,
	#a diagonal element is zero, it means rank<size of matrix.
	if not div_row:
	print "Check the circuit. Some variables have no dynamics."
	else:
	for c2 in range(matrix_e.columns):
	matrix_e.data[c1][c2]=matrix_e.data[c1][c2]/div_row
	for c2 in range(matrix_a.columns):
	matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/div_row
	for c2 in range(matrix_b.columns):
	matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/div_row
	#Row operations to remove row entries below the diagonal unity element.
	for c2 in range(c1+1, matrix_e.rows):
	mul_column=matrix_e.data[c2][c1]
	for c3 in range(matrix_e.columns):
	matrix_e.data[c2][c3]=matrix_e.data[c2][c3]-matrix_e.data[c1][c3]*mul_column
	for c3 in range(matrix_a.columns):
	matrix_a.data[c2][c3]=matrix_a.data[c2][c3]-matrix_a.data[c1][c3]*mul_column
	for c3 in range(matrix_b.columns):
	matrix_b.data[c2][c3]=matrix_b.data[c2][c3]-matrix_b.data[c1][c3]*mul_column
	return
	def mat_ode(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, runge_vectors):
	""" Solves the ODE using Runge Kutta 4th order method. """
	mat_ode_reduce(matrix_e, matrix_a, matrix_b)
	def runge_function(x_plus_dxdt, ddt_mat):
	""" Defines the function dx/dt=f(x) """
	#Calculate d/dt vector in reverse.
	for c1 in range(matrix_e.rows-1, -1, -1):
	if matrix_e.data[c1][c1]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	ddt_mat.data[c1][0]=matrix_b.data[c1][0]*input_vector
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	ddt_mat.data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	for c3 in range(matrix_a.columns):
	ddt_mat.data[c1][0]-=matrix_a.data[c1][c3]*x_plus_dxdt.data[c3][0]
	for c3 in range(c1+1, matrix_e.columns):
	ddt_mat.data[c1][0]-=matrix_e.data[c1][c3]*ddt_mat.data[c3][0]
	else:
	#If rank<size, set the (size-rank) solutions to zero.
	print "This result will be wrong!"
	ddt_mat.data[c1][0]=0.0
	return ddt_mat
	# Runge-Kutta 4th order method
	k1=runge_vectors[0]
	k1(runge_function(state_vector[0], runge_vectors[4])*dt) #k1=dt*(dx/dt)
	k2=runge_vectors[1]
	k2(runge_function(state_vector[0]+0.5*k1, runge_vectors[4])*dt) #k2=dt*d/dt(x+k1/2)
	k3=runge_vectors[2]
	k3(runge_function(state_vector[0]+0.5*k2, runge_vectors[4])*dt) #k3=dt*d/dt(x+k2/2)
	k4=runge_vectors[3]
	k4(runge_function(state_vector[0]+k3, runge_vectors[4])*dt) #k4=dt*d/dt(x+k3)
	for c1 in range(state_vector[1].rows):
	state_vector[1].data[c1][0]=(state_vector[0].data[c1][0]+(k1.data[c1][0]+2.0*k2.data[c1][0]+2.0*k3.data[c1][0]+k4.data[c1][0])*(1/6.0))
	R1=10.0
	R2=10.0
	R3=10.0
	L1=20.0e-3
	L2=20.0e-3
	L3=20.0e-3
	a=matrix.Matrix(2,2)
	b=matrix.Matrix(2)
	b.data=[[1.0],[0.0]]
	a.data=[[R1+R2, -R2], [-R2, R2+R3]]
	e=matrix.Matrix(2,2)
	e.data=[[L1+L2, -L2], [-L2, L2+L3]]
	curr_state_vec = matrix.Matrix(2)
	curr_state_vec.zeros(2)
	next_state_vec = matrix.Matrix(2)
	next_state_vec.zeros(2)
	dt=20.0e-6
	f=open("data_file.dat","w")
	input_volt=matrix.Matrix(2)
	ode_k1=matrix.Matrix(2)
	ode_k2=matrix.Matrix(2)
	ode_k3=matrix.Matrix(2)
	ode_k4=matrix.Matrix(2)
	ode_dbydt=matrix.Matrix(2)
	ode_solver=[ode_k1, ode_k2, ode_k3, ode_k4, ode_dbydt]
	for time_step in range(0,10000):
	t=time_step*dt
	u=100.0*math.sin(2*60*math.pi*t)
	mat_ode(e, a, b, [curr_state_vec, next_state_vec], u, 20.0e-6, ode_solver)
	curr_state_vec=next_state_vec
	f.write("%s \t %s \t %s \t %s \n" %(str(t), str(u), str(next_state_vec.data[0][0]), str(next_state_vec.data[1][0])))
	f.close()

view raw ode1.py hosted with ❤ by GitHub

It works with the basic resistor, inductor, voltage source:

 -------R1----L1---------------R3----L3------
|                             |                                 |
|                             |                                 |
|                             R2                               |
 V                            |                                 |
|                             L2                                |
|                              |                                 |
|                              |                                 |
|----------------------|-------------------------|

I hate drawing circuits. I draw so many of them while publishing papers, I can't get myself to draw for this blog. But I think this should give an idea of what I am simulating.

A few things to consider:

1. The code has been written with a consideration for simulation speed. For example, as far as possible, the arguments passed to functions are as references and are used as such. Copy contructors are not invoked and these matrices are changed in the function and reflect in the main program which calls the function.
The reason, if I were to make the functions black bloxes which accept inputs without modifying them and return outputs to be used in the larger algorithm, I would need to copy the matrices each time the function is called. This would mean, copy constructors and destructors called with each function call. A simulation would run for 1 second with a time step of 20 microseconds, which means hundreds of thousands of function calls with the same number of times the constuctors and destructors are called. This would give a normal desktop a severe bellyache.
So, by using references wherever possible, the simulation runs faster as there are minimum number of copy constructors.

2. The function call:
mat_ode(e, a, b, [curr_state_vec, next_state_vec], u, 20.0e-6, ode_solver)
Looks so neat but hides the fact that several vectors are grouped together and passed as a single object:
ode_solver=[ode_k1, ode_k2, ode_k3, ode_k4, ode_dbydt]
Here I must say that I like this feature of Python. You pass anything into a function and it is treated as an object. As long as you know what to do inside the function, you are OK. Moreover, here "ode_solver" is a mutable object - a list of lists (with lists)! So inside the function:
 k1=runge_vectors[0]
Will help me unpack the different components.

3. This means I need to take a look a the matrix.py class definition and try to make sure constructors are avoided as far as possible.

4. The ODE solves assuming that E is non-singular. But a look at the circuit will tell you that if L2=L3=0, E will be non-singular but will still be a valid circuit with one variable having a d/dt. So the ODE solver will have to check for singularity and solve spearately as a special case.


Anyway, this my guess will be the last post in 2012.

Matrix operations

So my first piece of code will be with respect to matrix operations. NumPy already has many in-built functions, but I still want to do this because I want to write my own code in the beginning rather than just use in-built components. So here goes matrix.py.

	#! /usr/bin/env python
	import sys
	class Matrix:
	""" For defining arrays of any given size.
	Contains overloaded operators for print, add, subtract and multiply.
	Does not contain = assignment operator! Use a(b) instead of a=b.
	Last updated - Dec 25, 2012 """
	def zeros(self, rows, columns=1):
	""" Creates a zero matrix
	Dimensions have to be passed as <Matrix>.zeros(rows, columns).
	This function also creates the actual array.
	So every method of the class/object should call this first before using the array."""
	self.rows=rows
	self.columns=columns
	self.data=[]
	for count1 in range(self.rows):
	row_vector=[]
	for count2 in range(self.columns):
	row_vector.append(0)
	self.data.append(row_vector)
	def __init__(self, rows=1, columns=1):
	"""Creates an object with a zero array of size rowsXcolumns.
	If only one dimension is specified specified, a column vector is assumed.
	If no dimensions are specified, an null matrix is created."""
	self.rows=rows
	self.columns=columns
	self.zeros(self.rows, self.columns)
	def __call__(self, matrix1):
	""" Replaces the assignment operator =. Use as c(b) instead of c=b.
	The = will only result in a reference to the object and not an actual copy.
	This method ensures it. """
	self.rows=matrix1.rows
	self.columns=matrix1.columns
	self.zeros(self.rows, self.columns)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	self.data[count1][count2]=matrix1.data[count1][count2]
	def identity(self, rows):
	""" Generates an identity matrix of rowsXrows. """
	self.rows=rows
	self.columns=rows
	self.zeros(self.rows, self.columns)
	for count in range(self.rows):
	self.data[count][count]=1.0
	def display(self):
	""" Explicit method for displaying an array. """
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	print self.data[count1][count2],
	print "\t",
	print "\n",
	def __repr__(self):
	""" Displays an array when print <Matrix> is used. """
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	print self.data[count1][count2],
	print "\t",
	print "\n",
	return " "
	def __add__(self, matrix1):
	""" Left addition of two matrices. """
	add_result=Matrix(self.rows, self.columns)
	try:
	if matrix1.rows: pass
	except AttributeError:
	print "*"*40
	print "Can't add a float or integer to a matrix directly."
	print "*"*40
	sys.exit("Logic Error")
	else:
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	print "Matrices not compatible for summation."
	sys.exit("Logic Error")
	else:
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	add_result.data[count1][count2]=self.data[count1][count2]+\
	matrix1.data[count1][count2]
	return add_result
	def __radd__(self, matrix1):
	""" Right addition of two matrices. """
	add_result=Matrix(self.rows, self.columns)
	try:
	if matrix1.rows: pass
	except AttributeError:
	print "*"*40
	print "Can't add a float or integer to a matrix directly."
	print "*"*40
	sys.exit("Logic Error")
	else:
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	print "Matrices not compatible for summation."
	else:
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	add_result.data[count1][count2]=self.data[count1][count2]+\
	matrix1.data[count1][count2]
	return add_result
	def __sub__(self, matrix1):
	""" Left subtraction of two matrices. """
	sub_result=Matrix(self.rows, self.columns)
	try:
	if matrix1.rows: pass
	except AttributeError:
	print "*"*40
	print "Can't subtract a float or integer to a matrix directly."
	print "*"*40
	sys.exit("Logic Error")
	else:
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	print "Matrices not compatible for subtraction."
	sys.exit("Logic Error")
	else:
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	sub_result.data[count1][count2]=self.data[count1][count2]-\
	matrix1.data[count1][count2]
	return sub_result
	def __rsub__(self, matrix1):
	""" Right subtraction of two mnatrices. """
	sub_result=Matrix(self.rows, self.columns)
	try:
	if matrix1.rows: pass
	except AttributeError:
	print "*"*40
	print "Can't subtract a float or integer to a matrix directly."
	print "*"*40
	sys.exit("Logic Error")
	else:
	if not (self.rows==matrix1.rows and self.columns==matrix1.columns):
	print "Matrices not compatible for subtraction."
	sys.exit("Logic Error")
	else:
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	sub_result.data[count1][count2]=matrix1.data[count1][count2]-\
	self.data[count1][count2]
	return sub_result
	def __mul__(self, matrix1):
	""" Left multiplication of two matrices.
	The multiplication is NOT done elementwise but in the conventional mathematical sense. """
	try:
	if matrix1.rows: pass
	except AttributeError:
	mul_result=Matrix(self.rows, self.columns)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	mul_result.data[count1][count2]=self.data[count1][count2]*matrix1
	else:
	mul_result=Matrix()
	if not (self.columns==matrix1.rows):
	print "Matrices not compatible for multiplication."
	else:
	mul_result.zeros(self.rows, matrix1.columns)
	for count1 in range(self.rows):
	for count2 in range(matrix1.columns):
	for count3 in range(self.columns):
	mul_result.data[count1][count2]+=self.data[count1][count3]*\
	matrix1.data[count3][count2]
	return mul_result
	def __rmul__(self, matrix1):
	""" Right multiplication of two matrices.
	The multiplication is NOT done elementwise but in the conventional mathematical sense. """
	try:
	if matrix1.rows: pass
	except AttributeError:
	mul_result=Matrix(self.rows, self.columns)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	mul_result.data[count1][count2]=self.data[count1][count2]*matrix1
	else:
	mul_result=Matrix()
	if not (self.rows==matrix1.columns):
	print "Matrices not compatible for multiplication."
	else:
	mul_result.zeros(matrix1.rows, self.columns)
	for count1 in range(matrix1.rows):
	for count2 in range(self.columns):
	for count3 in range(self.rows):
	mul_result.data[count1][count2]+=self.data[count1][count3]*\
	matrix1.data[count3][count2]
	return mul_result
	def transpose(self):
	""" Returns a matrix that is transpose of the calling object array.
	A new object is returned. Original object is untouched. """
	trans_result=Matrix(self.columns, self.rows)
	for count1 in range(self.rows):
	for count2 in range(self.columns):
	trans_result.data[count2][count1]=self.data[count1][count2]
	return trans_result

view raw matrix.py hosted with ❤ by GitHub

A few thoughts particularly with respect to my ex-darling C++.

1. There isn't the possibility of overloading the "=" assignment operator. So is I want "a" and "b" as arrays and I want to do:
a=b
This will only cause a reference "a" to point to the the object "b". A copy will not be made. So if I change "b" later, "a" changes too. Which is not what I want. So I need to define the __call__ method. This was I could do:
a(b)
But of course you would need to define "a" as a matrix before this is run.

2. Indentation is a requirement in Python. It improves the readability of code. True. But doing away with braces {}, or the begin/end statement blocks is probably going too far. Just a beginner's thought anyway.

3. Function/Operator overloading: In C++, you can overload a function any number of times with a different number of variables, different type of variables etc. So an overloaded function could accept a float or a string and two completely different functions can be written. In Python, everything is an object. So you can't specify float or Matrix before an argument in a function/method. You have to figure out what that is inside the functions. This I feel is unnecessary ambiguity.
I was trying to code the "*" oveload operator for matrices, and would like to have the possibility of multiplying two matrices or a matrix and a float/integer. But the only was to do that was inside the overloaded operator. A little googling told me that type checking was not a good idea. Probably because they want to have  the possibility fo redefining the types or maybe defining additional super-types in later releases of Python. So the only thing I could use was exception handling. The code works but this isn't really an exception. This is a normal varation to the manner in which the operator is used on matrices. So seems odd.

4. Two more functions may be written depending on need - inverse and upper triangularization or lower triangularization. This will be decided later. The idea would be solve equations of the order Ax=b or E*d/dt(x)=Ax+bu. Using inverse to solve equations is usually avoided as matrix even close to singular will cause the simulation to blow up. So more on this later.

The beginning

Well, not really the beginning. Started this two years ago but had to stop because I was too busy. Besides, I was in a company where a technical blog could have been a conflict of interest. Now back to a university setting, I can restart this extremely interesting activity.

I started coding for simulations in the year 2007. To some extent, the real reason was I was more or less fed up with the commercial circuit simulators. They were expensive and worse buggy. Another problem was that simulating complex circuits for long duration led to programs that ran for a long time. In a Windows setting that essentially means your computer is unusable during the simulation and sometimes need a restart after the simulations ends.

When I first started coding in C, the pain was incredible. The smallest blocks that were available in simulators had to be coded and then made into usable functions. But as the complexity of the circuits grew, I found that circuits with four or five converters would be simulated with never the fear of crashing. I could run six C programs simultaneously in independent terminals in Linux which could take hours and at the same time use my computer as if nothing was running in the background.

This led me to convert my C code to C++.  An immensely rewarding experience that finally ended with me completing my PhD. Anyway, since then I started learning Python. Python simplifies life to some extent with in-built functions. But most importantly, it is free and open source. With the rapid development taking place in Python, I thought this might be the right time to jump in.

One of the mistakes I made while coding in C/C++ was that I never documented my steps. At that time, I wasn't into blogging. Well now I am. So here goes. I am going to start coding in Python with the objective of building a circuit simulator specifically for power electronics. This blog is going to be free-flowing, so it sometimes may seem like to ramblings of a drunkard. Well I want it to be so.

So this begins this journey of a tech blog. I hope this lasts as long as I am in university as a post doctoral researcher which I hope will be for the next two years. So God save me and God save those who read this blog :)