Circuit Solver - V

A little progress. At least the ODE is solvable now. I used the strategy of row operations to get rid of the elements with low (or zero) time constants from as many rows as possible. This is the code (click on "view raw" below the code box to see it in a new window):

	# Getting rid of elements with small time constants
	# that appear multiple times in a circuit.
	for c1 in range(matrix_a.rows):
	for c2 in range(c1+1, matrix_a.columns):
	if matrix_a.data[c1][c2]:
	# Check is an off-diagonal element has a
	# small L/R ratio.
	if (matrix_e.data[c1][c2]/matrix_a.data[c1][c2])<0.1*dt:
	# Check if the diagonal element is row c1 has
	# smaller L/R than diagonal element in row c2
	# The idea is to arrange the equations as per
	# time constants. So manipulate the larger time
	# constants with respect to the smaller time constants
	# and ignore the ones that are below a threshold.
	if (matrix_e.data[c1][c1]/matrix_a.data[c1][c1])<(matrix_e.data[c2][c2]/matrix_a.data[c2][c2]):
	for c3 in range(matrix_a.columns):
	matrix_a.data[c2][c3]-=matrix_a.data[c1][c3]
	matrix_e.data[c2][c3]-=matrix_e.data[c1][c3]
	for c3 in range(matrix_b.columns):
	matrix_b.data[c2][c3]-=matrix_b.data[c1][c3]
	else:
	for c3 in range(matrix_a.columns):
	matrix_a.data[c1][c3]-=matrix_a.data[c2][c3]
	matrix_e.data[c1][c3]-=matrix_e.data[c2][c3]
	for c3 in range(matrix_b.columns):
	matrix_b.data[c1][c3]-=matrix_b.data[c2][c3]
	# Getting rid of the time constants that
	# are too small for the simulation time
	# step. Can generate a warning in a log report.
	for c1 in range(matrix_e.rows):
	if matrix_a.data[c1][c1]:
	if (matrix_e.data[c1][c1]/matrix_a.data[c1][c1])<0.1*dt:
	if matrix_e.data[c1][c1]:
	for c2 in range(matrix_a.columns):
	matrix_a.data[c1][c2]=matrix_a.data[c1][c2]/matrix_e.data[c1][c1]
	for c2 in range(matrix_b.columns):
	matrix_b.data[c1][c2]=matrix_b.data[c1][c2]/matrix_e.data[c1][c1]
	matrix_e.data[c1][c1]=0.0

view raw stiff_solver1.py hosted with ❤ by GitHub

So for a circuit:

The matrices are transformed as follows:

Original matrices are:
------------------------------------------------------------------------------------------------------------------
e=
0.0  0.0  0.0  0.0
0.0  0.2  0.0  0.1
0.0  0.0  0.0001  0.0001
0.0  0.1  0.0001  0.2001

a=
10000000.0  0.0  10000000.0  10000000.0
0.0  20.0  0.0  10.0
10000000.0  0.0  10000010.1  10000000.1
10000000.0  10.0  10000000.1  10000020.1

b=
1.0
0.0
0.0
0.0
------------------------------------------------------------------------------------------------------------------

These are manipulated to:
------------------------------------------------------------------------------------------------------------------
e=
0.0  0.0  0.0  0.0
0.0  0.2  0.0  0.1
0.0  0.0  0.0001  0.0001
0.0  0.1  0.0001  0.2001

a=
10000000.0  0.0  10000000.0  10000000.0
0.0  20.0  0.0  10.0
0.0  0.0  10.0999999996  0.0999999996275
0.0  10.0  0.0999999996275  20.0999999996

b=
1.0
0.0
-1.0
-1.0
------------------------------------------------------------------------------------------------------------------

With this the ODE solves with a time step of 10.0e-6 seconds. However, the stiff equation still is not solved satifsfactorily because of the entire row of large resistances. The way, loop current i1 would be calculated would be:

i1=(-10000000.0*i3-10000000.0*i4)/10000000.0

This would cause i1 to be of the same order as the other loop currents as this reduces to:

i1=-i3-i4

Which is wrong.

So this makes me wonder, why are the off-diagonal terms present in that case? Why not just do:

i1=u/10000000.0

Seems like that might be the solution. Because all I need is the sum total of all the resistance between the input signal in the loop to give me the approximate current. But this is then eventually ONLY the approximate current.