User reported issue: short circuit created with switch open

[oliver-phet]

A user wrote in with this bug and I was able to quickly reproduce:

Test device Dell Latitude E6530 Dell Chromebook P26T

Operating System Windows 7 Chrome OS

Browser Chrome

Problem description A certain configuration leads to an inaccurate short circuit when a switch is opened.

Steps to reproduce In intro mode create a circuit in this configuration: Positive terminal of battery ->wire ->switch ->light bulb ->wire ->eraser ->wire ->negative terminal of battery

With the switch closed, nothing happens (as expected). When the switch is opened the bulb lights up and the battery and eraser catch fire. This should not happen since no current should move as the circuit is open.

The same thing happens when a dollar bill is used.

If the current is investigated with the Ammeter it shows differing current throughout which also shouldn't happen.

Visuals

[oliver-phet]

In some quick testing on my end, it appears to possibly be related to the eraser.

[arouinfar]

I built a few different similar circuits (including the one below), but couldn't reproduce the issue.

The screenshot from the user is exhibiting the bunched-up current that was also reported in https://github.com/phetsims/circuit-construction-kit-common/issues/281 which makes me think it's somehow related.

Assigning to @samreid to take a look.

[oliver-phet]

Oh weird - I was immediately able to reproduce, but now when I try, it seems more elusive.

Here's a circuit I just built that displays the bug with no electron bunching:

[oliver-phet]

Actually, there does seem to be some bunching. Here's that same circuit I built:

You can see the current arrows and electrons move through the bulb, but then seem to disappear/stack up.

[oliver-phet]

And to be clear - there are no "0 length" wires or overlapped contacts in my circuit. And no component voltages/resistances were adjusted from default.

[arouinfar]

So strange! Thanks for investigating @oliver-phet. The electons are missing from the bulb in https://github.com/phetsims/circuit-construction-kit-common/issues/693#issuecomment-828817870 so there might be some bunching going on there too.

[samreid]

Some observations:

• I was able to reproduce the problem 3-4 times.
• I reproduced the problem with a dollar bill in place of the eraser (noting they have the same resistance of 1000000000)
• Then I went through a phase where, following the same procedure as above, I was not able to reproduce the problem about 5 times in a row.
• Next, I saw a different symptom, where an eraser in an open circuit flickered on fire, but only while dragging (other symptoms were not present).

I don't know what is causing this, but hypothesized that it is a numerical error with 1000000000 resistance. I tried lowering it by a factor of 10 and saw the same problem, then by a factor of 100 and saw the same symptom with lower current. But I identified a better way to model true insulators--we can omit them from the circuit solution as we do for open switches. I don't recall if we wanted erasers and dollars to be pure insulators or just very high resistors. But I'll try modeling them as pure insulators to see what happens. But I'm skeptical I will be able to test this accurately, since there are some times where I test 5 times in a row (with no code changes) and cannot see the problem.

• In order to make sure I could see if the proposal had a positive impact, I ran the experiment with no code changes 10 times. I was not able to see the error in any of the cases. I also tried higher and lower resistance is for the eraser, and I tried higher voltage is for the battery. The bug remained elusive.

I'll make the change that models eraser and dollar bill as pure insulators, but I'm going to need help figuring out this issue.

[samreid]

@arouinfar can you please review my notes & the commit, and advise?

[samreid]

@arouinfar pointed out that it reads 0 Ohms, that should be corrected.

[samreid]

We would like to revert the prior commit, since we would like the resistance to be non-infinite.

[samreid]

After that, we will ask @KatieWoe and QA to investigate this issue, see if it is possible to reliably reproduce.

@kathy-phet reminds @samreid to put it in the issues to test/verify column on the QA project board.

[KatieWoe]

This seems to happen more noticeably with higher voltages, and doesn't seem to need a 0 length wire, though it does make it easier.

[KatieWoe]

Clarification: The above was in published

[KatieWoe]

I can have trouble reproducing this with non-lightbulbs, but it still happens in master. Turning up voltage helps.

[samreid]

I reverted the change that omits eraser and dollar from the circuit matrix.

[samreid]

After some iteration, I was able to get to this error case. Note the current changes as you adjust the wire length:

[samreid]

After creating the above movie, I tried removing a few components to simplify the circuit. I was able to remove several parts and keep the buggy behavior, but then after removing something else (a wire?) it could no longer produce the problem at all, even when adding everything back.

UPDATE: In a later trial, I got a movie of it:

[samreid]

Here's a case where attaching the light bulb triggers the problem:

[samreid]

For this circuit, here are the equations:

Debugging circuit: 
resistors:node1 -> node2 @ 2.07037670465044e-10 Ohms,node2 -> node3 @ 10 Ohms,node4 -> node1 @ 0.000001 Ohms, 
batteries: node0 -> node4 @ 100000 Volts, 
currentSources: 

V0=0
1000000*V1+-1000000*V4+-4830038889.801163*V2+4830038889.801163*V1=0
4830038889.801163*V2+-4830038889.801163*V1+-0.1*V3+0.1*V2=0
0.1*V3+-0.1*V2=0
-I0_4+-1000000*V1+1000000*V4=0
-V0+V4=100000

A=
dim: 6x6
0 1 0 0 0 0 
0 0 4831038889.801163 -4830038889.801163 0 -1000000 
0 0 -4830038889.801163 4830038889.901163 -0.1 0 
0 0 0 -0.1 0.1 0 
-1 0 -1000000 0 0 1000000 
0 -1 0 0 0 1 

z=
dim: 6x1
0 
0 
0 
0 
0 
100000 

unknowns=
I0_4
V0
V1
V2
V3
V4

x=
dim: 6x1
0.02496337890625 
0 
99999.99999997504 
99999.99999997503 
99999.99999997504 
100000 

Solving this system by hand, I see that we should have obtained:

I0_4=0
V0=0
V1=V2=V3=V4=100000 

So this seems to be a numerical matrix solve error.

[samreid]

Using bignumber.js v9.0.1, I adapted LU Decomposition to be arbitrary-precision, then with this patch:

```js Index: js/LUDecomposition.js IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/LUDecomposition.js b/js/LUDecomposition.js --- a/js/LUDecomposition.js (revision c63a40f4ec3f0c86260c4fe18a7708079925afe5) +++ b/js/LUDecomposition.js (date 1628857320958) @@ -11,6 +11,8 @@ const ArrayType = window.Float64Array || Array; +BigNumber.config({ DECIMAL_PLACES: 1000 }) + class LUDecomposition { constructor( matrix ) { let i; @@ -20,7 +22,12 @@ this.matrix = matrix; // TODO: size! - this.LU = matrix.getArrayCopy(); + // this.LU = matrix.entries.map(entry=>new BigNumber(entry)) + this.LU = []; + matrix.entries.forEach(entry=>{ + this.LU.push(new BigNumber(entry)) + }); + const LU = this.LU; this.m = matrix.getRowDimension(); const m = this.m; @@ -31,7 +38,7 @@ this.piv[ i ] = i; } this.pivsign = 1; - const LUcolj = new ArrayType( m ); + const LUcolj = new Array(m); // Outer loop. @@ -39,7 +46,7 @@ // Make a copy of the j-th column to localize references. for ( i = 0; i < m; i++ ) { - LUcolj[ i ] = LU[ matrix.index( i, j ) ]; + LUcolj[ i ] = new BigNumber(LU[ matrix.index( i, j ) ]); } // Apply previous transformations. @@ -47,13 +54,16 @@ for ( i = 0; i < m; i++ ) { // Most of the time is spent in the following dot product. const kmax = Math.min( i, j ); - let s = 0.0; + let s =new BigNumber(0); for ( k = 0; k < kmax; k++ ) { const ik = matrix.index( i, k ); - s += LU[ ik ] * LUcolj[ k ]; + const a = new BigNumber(LU[ik]); + const b = LUcolj[k]; + s = s.plus(a.times(b)); + // s += LU[ ik ] * LUcolj[ k ]; } - LUcolj[ i ] -= s; + LUcolj[ i ] = LUcolj[ i ].minus(s) LU[ matrix.index( i, j ) ] = LUcolj[ i ]; } @@ -61,7 +71,7 @@ let p = j; for ( i = j + 1; i < m; i++ ) { - if ( Math.abs( LUcolj[ i ] ) > Math.abs( LUcolj[ p ] ) ) { + if ( LUcolj[ i ].abs() > LUcolj[ p ].abs() ) { p = i; } } @@ -83,7 +93,7 @@ if ( j < m && LU[ this.matrix.index( j, j ) ] !== 0.0 ) { for ( i = j + 1; i < m; i++ ) { - LU[ matrix.index( i, j ) ] /= LU[ matrix.index( j, j ) ]; + LU[ matrix.index( i, j ) ] = LU[ matrix.index( i, j ) ].dividedBy(LU[ matrix.index( j, j ) ]); } } } @@ -97,7 +107,7 @@ isNonsingular() { for ( let j = 0; j < this.n; j++ ) { const index = this.matrix.index( j, j ); - if ( this.LU[ index ] === 0 ) { + if ( this.LU[ index ].isZero() ) { return false; } } @@ -108,44 +118,44 @@ * @public * * @returns {Matrix} - */ - getL() { - const result = new Matrix( this.m, this.n ); - for ( let i = 0; i < this.m; i++ ) { - for ( let j = 0; j < this.n; j++ ) { - if ( i > j ) { - result.entries[ result.index( i, j ) ] = this.LU[ this.matrix.index( i, j ) ]; - } - else if ( i === j ) { - result.entries[ result.index( i, j ) ] = 1.0; - } - else { - result.entries[ result.index( i, j ) ] = 0.0; - } - } - } - return result; - } - - /** - * @public - * - * @returns {Matrix} - */ - getU() { - const result = new Matrix( this.n, this.n ); - for ( let i = 0; i < this.n; i++ ) { - for ( let j = 0; j < this.n; j++ ) { - if ( i <= j ) { - result.entries[ result.index( i, j ) ] = this.LU[ this.matrix.index( i, j ) ]; - } - else { - result.entries[ result.index( i, j ) ] = 0.0; - } - } - } - return result; - } + // */ + // getL() { + // const result = new Matrix( this.m, this.n ); + // for ( let i = 0; i < this.m; i++ ) { + // for ( let j = 0; j < this.n; j++ ) { + // if ( i > j ) { + // result.entries[ result.index( i, j ) ] = this.LU[ this.matrix.index( i, j ) ]; + // } + // else if ( i === j ) { + // result.entries[ result.index( i, j ) ] = 1.0; + // } + // else { + // result.entries[ result.index( i, j ) ] = 0.0; + // } + // } + // } + // return result; + // } + // + // /** + // * @public + // * + // * @returns {Matrix} + // */ + // getU() { + // const result = new Matrix( this.n, this.n ); + // for ( let i = 0; i < this.n; i++ ) { + // for ( let j = 0; j < this.n; j++ ) { + // if ( i <= j ) { + // result.entries[ result.index( i, j ) ] = this.LU[ this.matrix.index( i, j ) ]; + // } + // else { + // result.entries[ result.index( i, j ) ] = 0.0; + // } + // } + // } + // return result; + // } /** * @public @@ -178,16 +188,16 @@ * * @returns {number} */ - det() { - if ( this.m !== this.n ) { - throw new Error( 'Matrix must be square.' ); - } - let d = this.pivsign; - for ( let j = 0; j < this.n; j++ ) { - d *= this.LU[ this.matrix.index( j, j ) ]; - } - return d; - } + // det() { + // if ( this.m !== this.n ) { + // throw new Error( 'Matrix must be square.' ); + // } + // let d = this.pivsign; + // for ( let j = 0; j < this.n; j++ ) { + // d *= this.LU[ this.matrix.index( j, j ) ]; + // } + // return d; + // } /** * @public @@ -209,12 +219,14 @@ // Copy right hand side with pivoting const nx = matrix.getColumnDimension(); const Xmat = matrix.getArrayRowMatrix( this.piv, 0, nx - 1 ); + const entries = []; + Xmat.entries.forEach(e => entries.push(new BigNumber(e))); // Solve L*Y = B(piv,:) for ( k = 0; k < this.n; k++ ) { for ( i = k + 1; i < this.n; i++ ) { for ( j = 0; j < nx; j++ ) { - Xmat.entries[ Xmat.index( i, j ) ] -= Xmat.entries[ Xmat.index( k, j ) ] * this.LU[ this.matrix.index( i, k ) ]; + entries[ Xmat.index( i, j ) ] = entries[ Xmat.index( i, j ) ].minus(entries[ Xmat.index( k, j ) ].times(this.LU[ this.matrix.index( i, k ) ])); } } } @@ -222,14 +234,17 @@ // Solve U*X = Y; for ( k = this.n - 1; k >= 0; k-- ) { for ( j = 0; j < nx; j++ ) { - Xmat.entries[ Xmat.index( k, j ) ] /= this.LU[ this.matrix.index( k, k ) ]; + entries[ Xmat.index( k, j ) ] = entries[ Xmat.index( k, j ) ].dividedBy(this.LU[ this.matrix.index( k, k ) ]); } for ( i = 0; i < k; i++ ) { for ( j = 0; j < nx; j++ ) { - Xmat.entries[ Xmat.index( i, j ) ] -= Xmat.entries[ Xmat.index( k, j ) ] * this.LU[ this.matrix.index( i, k ) ]; + entries[ Xmat.index( i, j ) ] = entries[ Xmat.index( i, j ) ].minus(entries[ Xmat.index( k, j ) ].times(this.LU[ this.matrix.index( i, k ) ])); } } } + + const efloats = entries.map(e=>e.toNumber()); + Xmat.entries = efloats; return Xmat; } } ```

We can turn up the precision to get a more accurate result:

15 decimal places:

dim: 6x1
0.095781237 
0 
99999.99999990422 
99999.9999999042 
100000.00023704463 
100000 

20 decimal places

dim: 6x1
6.1854304e-7 
0 
100000 
100000 
100000.00000000368 
100000 

30 decimal places

dim: 6x1
-2.5093793e-17 
0 
100000 
100000 
100000 
100000 

100 decimal places

dim: 6x1
1.01516598e-86 
0 
100000 
100000 
100000 
100000 

1000 decimal places

dim: 6x1
0 
0 
100000 
100000 
100000 
100000 

[samreid]

By the way, here is the test harness I have been using:

const stringToMatrix = string => {
  const data = string.split( '\n' ).map( row => row.trim().split( ' ' ) );
  const m = data.length;
  const n = data[ 0 ].length;
  const matrix = new Matrix( m, n );
  for ( let i = 0; i < m; i++ ) {
    for ( let k = 0; k < n; k++ ) {
      matrix.set( i, k, parseFloat( data[ i ][ k ] ) );
    }
  }
  return matrix;
};
const matrix = stringToMatrix( `0 1 0 0 0 0 
0 0 4831038889.801163 -4830038889.801163 0 -1000000 
0 0 -4830038889.801163 4830038889.901163 -0.1 0 
0 0 0 -0.1 0.1 0 
-1 0 -1000000 0 0 1000000 
0 -1 0 0 0 1` );

console.log( matrix.toString() );

const b = stringToMatrix( `0 
0 
0 
0 
0 
100000 ` );
console.log( b.toString() );

const x = matrix.solve( b );
console.log( x.toString() );

[samreid]

@jonathanolson and I discussed the proposed approach, and agreed it seems reasonable to move forward with decimal.js in sherpa, and adding HighPrecisionLUDecomposition in dot, which would be a copy-pasted version of LUDecomposition (the parts needed for a solve), so both can run at full speed. We also tested and saw all matrices to be square, so we don't think we need to do QR (but if a non-square shows up, we will have to do QR). HighPrecisionLUDecomposition will use decimal.js, but only sims that need HighPrecisionLUDecomposition will have to include decimal.js since it is 40KB+ minified.

I am hopeful it will solve this corner case, and it remains to be seen if it may solve other corner cases as well.

[samreid]

I have a patch with decimal.js working nicely, but in testing before I committed, I observed that even with 200 decimal places, you can easily create a circuit that has a bad solution:

Not sure if I should commit the proposal as a partial solution, or whether we need to find another linear algebra solution.

[samreid]

I pushed what @jonathanolson and I discussed even though I can see it still fails on some easy cases. Maybe there is still a bug or imprecise part in LUDecompositionDecimal, or maybe we need a new algorithm.

[samreid]

Here's a case that's still failing:

Debugging circuit: resistors:node1 -> node2 @ 2.0915055657130383e-10 Ohms,node2 -> node3 @ 10 Ohms,node4 -> node1 @ 0.000001 Ohms, batteries: node0 -> node4 @ 100000 Volts, currentSources: 
    equations:
    V0=0
1000000*V1+-1000000*V4+-4781244747.2932205*V2+4781244747.2932205*V1=0
4781244747.2932205*V2+-4781244747.2932205*V1+-0.1*V3+0.1*V2=0
0.1*V3+-0.1*V2=0
-I0_4+-1000000*V1+1000000*V4=0
-V0+V4=100000

    A=
dim: 6x6
0 1 0 0 0 0 
0 0 4782244747.2932205 -4781244747.2932205 0 -1000000 
0 0 -4781244747.2932205 4781244747.393221 -0.1 0 
0 0 0 -0.1 0.1 0 
-1 0 -1000000 0 0 1000000 
0 -1 0 0 0 1 

    z=
dim: 6x1
0 
0 
0 
0 
0 
100000 

    unknowns=
I0_4
V0
V1
V2
V3
V4

    x=
dim: 6x1
0.049999999999974995 
0 
99999.99999995 
99999.99999994999 
99999.99999994999 
100000 

Solving it by hand, it should still be 0.0. Either LU is not suited for this problem, or we have a bug in our LU.

[samreid]

I used math.js and lusolve, but it had the same bad error we started with:

```diff Index: js/model/ModifiedNodalAnalysisCircuit.js IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/ModifiedNodalAnalysisCircuit.js b/js/model/ModifiedNodalAnalysisCircuit.js --- a/js/model/ModifiedNodalAnalysisCircuit.js (revision b2a91530f36d82cefb497d36afbb9ba66c7b12cb) +++ b/js/model/ModifiedNodalAnalysisCircuit.js (date 1628888249706) @@ -330,43 +330,68 @@ for ( let i = 0; i < equations.length; i++ ) { equations[ i ].stamp( i, A, z, getIndex ); } - phet.log && phet.log( A.m, A.n ); - phet.log && phet.log( `Debugging circuit: ${this.toString()}` ); - phet.log && phet.log( equations.join( '\n' ) ); - phet.log && phet.log( `A=\n${A.toString()}` ); - phet.log && phet.log( `z=\n${z.toString()}` ); - phet.log && phet.log( `unknowns=\n${unknowns.map( u => u.toTermName() ).join( '\n' )}` ); // solve the linear matrix system for the unknowns let x; - try { - assert && assert( A.m === A.n, `the matrix should be square, instead it was ${A.m} x ${A.n}` ); - x = new LUDecompositionDecimal( A, LUDecimal ).solve( z, LUDecimal ); + // try { + assert && assert( A.m === A.n, `the matrix should be square, instead it was ${A.m} x ${A.n}` ); + x = new LUDecompositionDecimal( A, LUDecimal ).solve( z, LUDecimal ); - // const x1 = A.inverse().times( z ); - // console.log( 'x' ) - // console.log( x.toString() ); - // console.log( 'x1' ) - // console.log( x1.toString() ); - // - // const a1 = new LUDecompositionDecimal( A, LUDecimal ).solve( Matrix.identity( A.m, A.m ), LUDecimal ); - // const x2 = a1.times( z ); - // console.log( 'x2' ); - // console.log( x2.toString() ); - } - catch( e ) { + // const x1 = A.inverse().times( z ); + // console.log( 'x' ) + // console.log( x.toString() ); + // console.log( 'x1' ) + // console.log( x1.toString() ); + // + // const a1 = new LUDecompositionDecimal( A, LUDecimal ).solve( Matrix.identity( A.m, A.m ), LUDecimal ); + // const x2 = a1.times( z ); + // console.log( 'x2' ); + // console.log( x2.toString() ); - // Sometimes a fuzz test gives a deficient matrix rank. It is a rare error and I haven't got one in the - // debugger yet to understand the cause. Catch it and provide a solution of zeroes of the correct dimension - // See https://github.com/phetsims/circuit-construction-kit-dc/issues/113 - x = new Matrix( A.n, 1 ); + console.log( A.toString() ); + const m = []; - // console.log( 'Rank deficient matrix!' ) + for ( let i = 0; i < A.m; i++ ) { + const e = A.getMatrix( i, i, 0, A.n - 1 ).entries; + const x = []; + e.forEach( element => { + x.push( element ); + } ); + m.push( x ); } + const zz = []; + z.entries.forEach( zzz => { + zz.push( zzz ); + } ); + const x3 = math.lusolve( m, zz); + console.log( 'x3' ); + console.log( x3 ); + // } + // catch( e ) { + // + // // Sometimes a fuzz test gives a deficient matrix rank. It is a rare error and I haven't got one in the + // // debugger yet to understand the cause. Catch it and provide a solution of zeroes of the correct dimension + // // See https://github.com/phetsims/circuit-construction-kit-dc/issues/113 + // x = new Matrix( A.n, 1 ); + // + // // console.log( 'Rank deficient matrix!' ) + // } // The matrix should be square since it is an exact analytical solution, see https://github.com/phetsims/circuit-construction-kit-dc/issues/96 assert && assert( A.m === A.n, 'Matrix should be square' ); - phet.log && phet.log( `x=\n${x.toString()}` ); + + console.log( `Debugging circuit: ${this.toString()} + equations: + ${equations.join( '\n' )} + + A=\n${A.toString()} + + z=\n${z.toString()} + + unknowns=\n${unknowns.map( u => u.toTermName() ).join( '\n' )} + + x=\n${x.toString()} + ` ); const voltageMap = {}; for ( let i = 0; i < unknownVoltages.length; i++ ) { ```

Matlab gives the same type of answer:

[samreid]

Here's the matrix problem:

a = 1000000
b = 4781244747.2932205
V0 = 0
a*V1 -a*V4 -b*V2 + b*V1 = 0
b*V2 -b*V1 -0.1*V3 + 0.1*V2 = 0
0.1*V3 -0.1*V2 = 0
-I0_4 -a*V1 + a*V4 = 0
-V0 + V4 = 100000

solving through substitution (by hand) gives I0_4 = 0, but all the numerical methods put it around +/-0.01 or so.

[samreid]

The other part that is unsolved is that the circuit mostly has the correct solution, but as the wire is resized it goes through spots of bad solution:

[samreid]

Anonymizing the values, here are snapshots of 2 pixels apart (bad on the left, good on the right):

The only difference is in the wire resistance. The solver for both is the LUDecompositionDecimal with 200 significant digits.

Patch with debugging output:

```diff Index: js/model/ModifiedNodalAnalysisCircuit.js IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/js/model/ModifiedNodalAnalysisCircuit.js b/js/model/ModifiedNodalAnalysisCircuit.js --- a/js/model/ModifiedNodalAnalysisCircuit.js (revision b2a91530f36d82cefb497d36afbb9ba66c7b12cb) +++ b/js/model/ModifiedNodalAnalysisCircuit.js (date 1628893913593) @@ -324,49 +324,83 @@ return index; }; + debugger; // Prepare the A and z matrices for the linear system Ax=z const A = new Matrix( equations.length, this.getNumVars() ); const z = new Matrix( equations.length, 1 ); + + console.log( 'starting A:' + A.toString() ); for ( let i = 0; i < equations.length; i++ ) { + console.log( 'stamping row ' + i ); equations[ i ].stamp( i, A, z, getIndex ); + // equations[ i ] + + console.log( 'stamped: ' + equations[ i ].toString() ); + console.log( 'A:' + A.toString() ); } - phet.log && phet.log( A.m, A.n ); - phet.log && phet.log( `Debugging circuit: ${this.toString()}` ); - phet.log && phet.log( equations.join( '\n' ) ); - phet.log && phet.log( `A=\n${A.toString()}` ); - phet.log && phet.log( `z=\n${z.toString()}` ); - phet.log && phet.log( `unknowns=\n${unknowns.map( u => u.toTermName() ).join( '\n' )}` ); // solve the linear matrix system for the unknowns let x; - try { - assert && assert( A.m === A.n, `the matrix should be square, instead it was ${A.m} x ${A.n}` ); - x = new LUDecompositionDecimal( A, LUDecimal ).solve( z, LUDecimal ); + // try { + assert && assert( A.m === A.n, `the matrix should be square, instead it was ${A.m} x ${A.n}` ); + debugger; + x = new LUDecompositionDecimal( A, LUDecimal ).solve( z, LUDecimal ); - // const x1 = A.inverse().times( z ); - // console.log( 'x' ) - // console.log( x.toString() ); - // console.log( 'x1' ) - // console.log( x1.toString() ); - // - // const a1 = new LUDecompositionDecimal( A, LUDecimal ).solve( Matrix.identity( A.m, A.m ), LUDecimal ); - // const x2 = a1.times( z ); - // console.log( 'x2' ); - // console.log( x2.toString() ); - } - catch( e ) { + // const x1 = A.inverse().times( z ); + // console.log( 'x' ) + // console.log( x.toString() ); + // console.log( 'x1' ) + // console.log( x1.toString() ); + // + // const a1 = new LUDecompositionDecimal( A, LUDecimal ).solve( Matrix.identity( A.m, A.m ), LUDecimal ); + // const x2 = a1.times( z ); + // console.log( 'x2' ); + // console.log( x2.toString() ); - // Sometimes a fuzz test gives a deficient matrix rank. It is a rare error and I haven't got one in the - // debugger yet to understand the cause. Catch it and provide a solution of zeroes of the correct dimension - // See https://github.com/phetsims/circuit-construction-kit-dc/issues/113 - x = new Matrix( A.n, 1 ); - - // console.log( 'Rank deficient matrix!' ) - } + // console.log( A.toString() ); + // const m = []; + // + // for ( let i = 0; i < A.m; i++ ) { + // const e = A.getMatrix( i, i, 0, A.n - 1 ).entries; + // const x = []; + // e.forEach( element => { + // x.push( element ); + // } ); + // m.push( x ); + // } + // const zz = []; + // z.entries.forEach( zzz => { + // zz.push( zzz ); + // } ); + // const x3 = math.lusolve( m, zz); + // console.log( 'x3' ); + // console.log( x3 ); + // } + // catch( e ) { + // + // // Sometimes a fuzz test gives a deficient matrix rank. It is a rare error and I haven't got one in the + // // debugger yet to understand the cause. Catch it and provide a solution of zeroes of the correct dimension + // // See https://github.com/phetsims/circuit-construction-kit-dc/issues/113 + // x = new Matrix( A.n, 1 ); + // + // // console.log( 'Rank deficient matrix!' ) + // } // The matrix should be square since it is an exact analytical solution, see https://github.com/phetsims/circuit-construction-kit-dc/issues/96 assert && assert( A.m === A.n, 'Matrix should be square' ); - phet.log && phet.log( `x=\n${x.toString()}` ); + + console.log( `Debugging circuit: ${this.toString()} + equations: + ${equations.join( '\n' )} + + A=\n${A.toString()} + + z=\n${z.toString()} + + unknowns=\n${unknowns.map( u => u.toTermName() ).join( '\n' )} + + x=\n${x.toString()} + ` ); const voltageMap = {}; for ( let i = 0; i < unknownVoltages.length; i++ ) { @@ -535,10 +569,10 @@ * @param {number} row - the index of the row for this equation * @param {Matrix} a - the matrix of coefficients in Ax=z * @param {Matrix} z - the matrix on the right hand side in Ax=z - * @param {function} getIndex - (UnknownCurrent|UnknownVoltage) => number + * @param {function} getColumn - (UnknownCurrent|UnknownVoltage) => number * @public */ - stamp( row, a, z, getIndex ) { + stamp( row, a, z, getColumn ) { // Set the equation's value into the solution matrix z.set( row, 0, this.value ); @@ -546,9 +580,9 @@ // For each term, augment the coefficient matrix for ( let i = 0; i < this.terms.length; i++ ) { const term = this.terms[ i ]; - const index = getIndex( term.variable ); - a.set( row, index, term.coefficient + a.get( row, index ) ); - assert && assert( !isNaN( term.coefficient + a.get( row, index ) ), 'term should be a number' ); + const column = getColumn( term.variable ); + a.set( row, column, term.coefficient + a.get( row, column ) ); + assert && assert( !isNaN( term.coefficient + a.get( row, column ) ), 'term should be a number' ); } } ```

[samreid]

@jonathanolson and I met about this today and confirmed the matrices are ill-conditioned. We used Mathematica and singularValueDecomposition[2] // MatrixForm and saw that it was typically 1E10 (high value) vs 0.01 (low value). @jonathanolson identified a spot where an addition is done while creating the matrix, where there is an opportunity for floating point error. We hypothesize that using arbitrary precision for that addition may stabilize the solution (even though the matrix remains ill-conditioned). I'll test that next.

[samreid]

I added the higher precision for term addition, and it seems to stabilize the problem. Even though the matrix is ill-conditioned, we arrive at a stable, more accurate result. This is great news! Testing with 200 points of precision in CCK-AC yields performance problems, so I think the next step would be to tune the precision/performance balance.

[samreid]

I defaulted the precision to 32 places, and added a query parameter "?precision" in case we want to tune it further.

[samreid]

I think this will be ready for testing in the next round. On hold until other issues are fixed.

[samreid]

While working on, https://github.com/phetsims/circuit-construction-kit-common/issues/678, I observed this:

[samreid]

I requested QA testing for this issue. Please try creating a variety of circuits like those in the preceding comments to see if the electrons move in an open circuit.

[stemilymill]

I found some instances on master (most recent built version)

[KatieWoe]

What was the day of the build?

[samreid]

The current timestamp of the cck ac build on phettest is 2021-08-24

[stemilymill]

I built a version today to continue testing, and found similar issues (all from same circuit, no reset in between)

[samreid]

Is there a set of steps that can be used to reliably reproduce this problem? Also, can you try with different values for ?precision=... to see if that reduces the problem? The default precision is 32, but you could try increasing it.

[stemilymill]

It happens very occasionally, so I haven't figured out a consistent way to reproduce it. I will try changing the values for precision query.

When the bug happens, it happens right after the particular circuit is built, there doesn't seem to be a need for repeated stretching of wires for instance.

[samreid]

The slow solver is causing problems elsewhere, I wonder if we should go back to the floating point numerical solution and instead increase the wire resistance. The present value of ?wireResistivity is 1E-12 which gives 1E-10 ohms for a long wire. How was this value determined and can we turn it up? For instance, ?wireResistivity=1E-8 gives 0.00009 Ohms but has better behavior under the circumstances above. We can also experiment with turning down ?precision, which is currently at a default value of ?precision=32. Keep in mind that ?precision=16 is equal to machine precision, but takes much longer to compute.

We could also do more complicated things, like running the solver with precision=16 and precision=20, comparing answers and if they are different, going to precision=32. If they are the same, we already have the answer. But that means running the solver more times per frame.

[kathy-phet]

The 1E-12 was determined specifically to solve some of the issues around parallel batteries in DC circuits. Can you find the related issue there, and link it here. We worked hard to hone the behavior there. I think we could change it by a factor of 10, but probably not more than that.

[samreid]

We changed the wire resistivity to 1E-12 in the commit following this comment: https://github.com/phetsims/circuit-construction-kit-common/issues/676#issuecomment-811012072

[samreid]

In today's design meeting, I described the solution that runs the high precision solver 10 times before turning to the high performance solver. Also, we saw a dangling circuit element that led to buggy behavior. I described the ?precision=... query parameter to QA so they can see if turning up the precision can avoid errors like that.

[Nancy-Salpepi]

Using version 1.0.0-dev.20 on MacBook Air + chrome, while in the process of creating a circuit, the electrons moved, the power source went on fire and the bulb lit (even though it was still an incomplete circuit).

Steps to reproduce the above:

1. Place bulb in play area followed by an AC voltage source
2. Connect a wire to the bulb and then to the AC voltage
3. Connect a wire on the right side of the bulb
4. Connect a wire on the left side of the bulb (where the first wire was placed)
5. Wiggle the circuit around

[samreid]

Thanks for the great instructions @Nancy-Salpepi, I can easily produce this problem. I tried with ?precision=1000 and still see the problem. I tried ?capacitorResistance=1E-4 and still see the problem. Maybe the solution is to eliminate "dangling" circuit elements from the matrix. If a circuit element is not a member of a loop, then it doesn't contribute equations to the matrix, and doesn't get values from the matrix.

[samreid]

In the commit, I made it so that only circuit elements in a loop with a voltage source contribute to the circuit solution. Any circuit element not in a loop with a voltage source will have its current zeroed out. This was a complicated depth-first-search, and has not been well-vetted. Open switches are not traversed. Capacitors and inductors contribute a voltage if their companion model provides a voltage.

To test this issue, we should check that the bugs above no longer occur, but also that we still get the expected results in other scenarios. Please test using version: https://phet-dev.colorado.edu/html/circuit-construction-kit-ac/1.0.0-dev.23/phet/circuit-construction-kit-ac_en_phet.html

I noticed that disconnecting a charged capacitor makes it look uncharged, but it still seems charged. I'll open a side issue for that one.

[Nancy-Salpepi]

@sam reid @stemilymill and I tried several several times to replicate the issues seen in https://github.com/phetsims/circuit-construction-kit-common/issues/693#issuecomment-906668899 and https://github.com/phetsims/circuit-construction-kit-common/issues/693#issuecomment-914612755 but could not.

Can you elaborate on how you would like me to test: "but also that we still get the expected results in other scenarios?"

[samreid]

Can you elaborate on how you would like me to test: "but also that we still get the expected results in other scenarios?"

Basically, testing other known scenarios, like:

• 2 resistors in series should add
• 2 resistors in parallel should have a harmonic sum like https://www.electronics-tutorials.ws/resistor/res_4.html
• battery resistor capacitor should have exponential decay with a time constant RC https://www.electronics-tutorials.ws/rc/rc_1.html
• switch behaves as expected
• etc