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.