Initializing new stiff equations

Another problem that arises when there are a large number of nonlinear devices is that nonstiff loops become stiff loops but their current has still not become negligible because the diode has turned off at the previous iteration. For example, current through a diode is -0.1A causing the diode to turn off and it's resistance to change to megaohms. When this stiff loops interacts with other stiff loops and the currents of those stiff loops are calculated, all other stiff loops end up having currents of close magnitudes of 0.1 A. Result: entire simulation is disrupted.

Two things need to be done:
1. Need to identify which loops have just become stiff.
2. Calculate those stiff loops before running the ode solver.

The code to identify the newly formed stiff loops is:

	# Check if the system is stiff and recalculate if it is.
	# The stiff loops are isolated to the minimum number of loops
	# So that the system dynamics are captured as far as possible.
	remove_stiffness(sys_mat_e, sys_mat_a, sys_mat_b, dt, branch_params, stiff_ratio, system_loop_map, source_list, [curr_state_vec, next_state_vec], ode_var, loop_stiff_info)
	# This list if of those loops that were nonstiff
	# and have become stiff in the previous calculation
	loop_event=[]
	for c1 in range(len(system_loop_map)):
	is_loop_stiff="no"
	loop_event.append("no")
	for c2 in range(len(system_loop_map[c1])):
	if system_loop_map[c1][c2]=="stiff_forward" or system_loop_map[c1][c2]=="stiff_reverse":
	is_loop_stiff="yes"
	if loop_stiff_info[c1]=="no" and is_loop_stiff=="yes":
	loop_event[c1]="yes"
	else:
	loop_event[c1]="no"
	# This is just a list of which loops are stiff and which are nonstiff
	if is_loop_stiff=="yes":
	loop_stiff_info[c1]="yes"
	else:
	loop_stiff_info[c1]="no"
	# Recalculate the newly formed stiff loops.
	mat_ode_stiff(sys_mat_e, sys_mat_a, sys_mat_b, [curr_state_vec, next_state_vec], sys_mat_u, t-t_ode_prev, ode_var, system_loop_map, node_list, loop_stiff_info, branch_events, loop_event)

view raw init_stiff_loops.py hosted with ❤ by GitHub

And running the stiff loop ode:

	def stiff_equation_solver(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, ode_vectors, sys_loop_map, list_of_nodes, loop_info_stiff, branch_info_event, ode_row):
	""" Solving stiff loops that are algebraic equations only """
	if not (matrix_a.data[ode_row][ode_row]):
	#The variable has no dynamics and can't even be calculated statically!
	# May be due to a redundant loop.
	ode_vectors[4].data[ode_row][0]=0.0
	state_vector[0].data[ode_row][0]=0.0
	state_vector[1].data[ode_row][0]=0.0
	else:
	ode_vectors[4].data[ode_row][0]=0.0
	state_vector[0].data[ode_row][0]=0.0
	state_vector[1].data[ode_row][0]=0.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[ode_row][0]+=matrix_b.data[ode_row][0]*input_vector
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][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[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	for c2 in range(matrix_a.columns):
	if not (ode_row==c2):
	state_vector[0].data[ode_row][0]-=matrix_a.data[ode_row][c2]*state_vector[1].data[c2][0]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	for c2 in range(matrix_e.columns):
	state_vector[0].data[ode_row][0]-=matrix_e.data[ode_row][c2]*ode_vectors[4].data[c2][0]/dt
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	state_vector[0].data[ode_row][0]=state_vector[0].data[ode_row][0]/matrix_a.data[ode_row][ode_row]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	return
	def mat_ode_stiff(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, ode_vectors, sys_loop_map, list_of_nodes, loop_info_stiff, branch_info_event, loop_events):
	""" Solves thoses loops that are stiff """
	for c1 in range(len(loop_info_stiff)-1, -1, -1):
	if loop_events[c1]=="yes":
	stiff_equation_solver(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, ode_vectors, sys_loop_map, list_of_nodes, loop_info_stiff, branch_info_event, c1)
	for c1 in range(len(loop_info_stiff)-1, -1, -1):
	if loop_info_stiff[c1]=="yes":
	stiff_equation_solver(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, ode_vectors, sys_loop_map, list_of_nodes, loop_info_stiff, branch_info_event, c1)
	return

view raw stiff_loop_solver.py hosted with ❤ by GitHub

One other thing still needs to be done. What if the newly formed stiff loops themselves are interacting with each other through a stiff branch? Now even this ode solver would fail because they would disrupt each other. The only thing to do would be do perform an upper triangularization of the newly formed stiff loops with respect to the newly formed stiff branches and solve them backwards. Will do this later.