## Sunday, March 10, 2013

### Stiff systems - II

Stupid mistake. The algorithm wasn't totally wrong. Still is wrong but not completely wrong. The problem was in the loop interaction.

Take this sample circuit:

The large resistors are Resistor_V2 (at 4D and of 500000.0 ohm) and Resistor_V1 (at 12 J and of 10000000.0 ohm).

The matrices turned out to be:
-----------------------------------------------------------------------------------------------------------
500000.0  500000.0  0.0  500000.0  500000.0
500000.0  500010.1  0.0  -500000.1  500000.1
0.0  0.0  10000010.0  0.0  -10000000.0
500000.0  -500000.1  0.0  500020.1  -499990.1
500000.0  500000.1  -10000000.0  -499990.1  10500010.1

0.0  0.0  0.0  0.0  0.0
0.0  0.0001  0.0  -0.0001  0.0001
0.0  0.0  0.1  0.0  0.0
0.0  -0.0001  0.0  0.2001  0.0999
0.0  0.0001  0.0  0.0999  0.1001

1.0
0.0
0.0
0.0
0.0
-----------------------------------------------------------------------------------------------------------

For loops defined as:
-----------------------------------------------------------------------------------------------------------
1D 2D 3D 4D 5D 6D 6C 6B 6A 5A 4A 3A 2A 1A 1B 1C 1D

1H 1G 1F 1E 1D 2D 3D 4D 5D 6D 6E 6F 6G 6H 5H 4H 3H 2H 1H

9L 10L 11L 12L 12K 12J 12I 12H 11H 10H 9H 9I 9J 9K 9L

1H 1G 1F 1E 1D 2D 3D 4D 5D 6D 6E 6F 6G 6H 6I 6J 6K 6L 5L 4L 3L 2L 1L 1K 1J 1I 1H

1H 1G 1F 1E 1D 2D 3D 4D 5D 6D 6E 6F 6G 6H 7H 8H 9H 10H 11H 12H 12I 12J 12K 12L 11L 10L 9L 8L 7L 6L 5L 4L 3L 2L 1L 1K 1J 1I 1H
-----------------------------------------------------------------------------------------------------------

Which is wrong. Because, the interaction between loop 2 and loop 4 is -500000.1 in the matrix A but should be positive from the sense of the loops above. The problem is in the code that defines the "forward" and "reverse" directions.

This is the new code (click on "view raw" to see the code in a new window):

The problem is with the judgement of length of branch with respect to the start and end nodes. The length of the branch will the number of elements in the branch. Which will be (end_node_st_node+1). This is exact block where there was the error:
-----------------------------------------------------------------------------------------------------------
if (abs(loop1_end_node-loop1_st_node)!= (len(sys_loop_off_diag[c3])-1)):
if (len(loop_branches[c1])-1-loop1_st_node==(len(sys_loop_off_diag[c3])-1)):
loop1_end_node=len(loop_branches[c1])-1
if (len(loop_branches[c1])-1-loop1_end_node==(len(sys_loop_off_diag[c3])-1)):
loop1_st_node=len(loop_branches[c1])-1
if (abs(loop2_end_node-loop2_st_node)!= (len(sys_loop_off_diag[c3])-1)):
if (len(loop_branches[c2])-1-loop2_st_node==(len(sys_loop_off_diag[c3])-1)):
loop2_end_node=len(loop_branches[c2])-1
if (len(loop_branches[c2])-1-loop2_end_node==(len(sys_loop_off_diag[c3])-1)):
loop2_st_node=len(loop_branches[c2])-1
-----------------------------------------------------------------------------------------------------------

So, the row manipulations to remove the stiff elements from the matrix produce the following matrices:
-----------------------------------------------------------------------------------------------------------
a=
500000.0  0.0  0.0  0.0  0.0
0.0  10.1  0.0  0.0999999999767  0.0999999999767
0.0  0.0  10500010.0  0.0  0.0
0.0  0.0999999999767  0.0  20.1  10.1
0.0  0.0999999999767  10.0  10.1  10.0999999996

e=
0.0  0.0  0.0  0.0  0.0
0.0  0.0001  0.0  0.0001  0.0001
0.0  0.0  0.0  0.0  0.0
0.0  0.0001  0.0  0.2001  0.1001
0.0  0.0001  0.1  0.1001  0.1001

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

The isolation has taken place. Because the two large resistors have been isolated to two separate rows. But the loop resistance between the 10000000.0 ohm resistor and the source is wrong. It should be 10000010.1 instead of 10500010.0.

Well, it doesn't blow up like before. So an improvement.

Releasing this as a new patch anyway (v0.1.2). Check on the SourceForge site soon.