Kirchhoff’s circuit laws and Z3 SMT-solver

The circuit I've created on falstad.com:

Click here to open it on their website and run: http://tinyurl.com/y8raoud3.

The problem: find all 3 current values in 2 loops. This is usually solved by solving a system of linear equations.

Overkill, but Z3 SMT-solver can be used here as well, since it can solve linear equations as well, over real numbers:

from z3 import *

i1, i2, i3 = Reals ("i1 i2 i3")

R1=2000
R2=5000
R3=1000

V1=5 # left
V2=15 # right

s=Solver()

s.add(i3 == i1+i2)

s.add(V1 == R1*i1 + R3*i3)
s.add(V2 == R2*i2 + R3*i3)

print s.check()
m=s.model()
print m
print m[i1].as_decimal(6)
print m[i2].as_decimal(6)
print m[i3].as_decimal(6)


And the result:

sat
[i3 = 11/3400, i1 = 3/3400, i2 = 1/425]
0.000882?
0.002352?
0.003235?

Same as on falstad.com online simulator.

Z3 represents real numbers as fractions, then we convert them to numerical form...

Further work: take a circuit as a graph and build a system of equations.

Gaussian elimination

SMT-solver is overkill, these linear equations can be solved using simple and well-known Gaussian elimination.

First, we rewrite the system of equation:

   i1 +    i2 -    i3 == 0
R1*i1 +         R3*i3 == V1
        R2*i2 + R3*i3 == V2

Or in matrix form:

[ 1,    1,    -1    | 0  ]
[ 2000, 0,    1000  | 5  ]
[ 0,    5000, 1000  | 15 ]

I can solve it using Wolfram Mathematica, using RowReduce:

In[1]:= RowReduce[{{1, 1, -1, 0}, {2000, 0, 1000, 5}, {0, 5000, 1000, 15}}]
Out[1]= {{1,0,0,3/3400},{0,1,0,1/425},{0,0,1,11/3400}}

In[2]:= 3/3400//N
Out[2]= 0.000882353

In[3]:= 1/425//N
Out[3]= 0.00235294

In[4]:= 11/3400//N
Out[4]= 0.00323529

This is the same result: i1, i2 and i3 in numerical form.

ReduceRow output is:

[ 1,0,0 | 3/3400  ]
[ 0,1,0 | 1/425   ]
[ 0,0,1 | 11/3400 ]

... back to expressions, this is:

1*i1 + 0*i2 + 0*i3 = 3/3400
0*i1 + 1*i2 + 0*i3 = 1/425
0*i1 + 0*i2 + 1*i3 = 11/3400

In other words, this is just what i1/i2/i3 are.

Now something down-to-earth, C example I've copypasted from Rosetta Code, working with no additional libraries, etc:

// copypasted from https://rosettacode.org/wiki/Gaussian_elimination#C

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
 
#define mat_elem(a, y, x, n) (a + ((y) * (n) + (x)))
 
void swap_row(double *a, double *b, int r1, int r2, int n)
{
	double tmp, *p1, *p2;
	int i;
 
	if (r1 == r2) return;
	for (i = 0; i < n; i++) {
		p1 = mat_elem(a, r1, i, n);
		p2 = mat_elem(a, r2, i, n);
		tmp = *p1, *p1 = *p2, *p2 = tmp;
	}
	tmp = b[r1], b[r1] = b[r2], b[r2] = tmp;
}
 
void gauss_eliminate(double *a, double *b, double *x, int n)
{
#define A(y, x) (*mat_elem(a, y, x, n))
	int i, j, col, row, max_row,dia;
	double max, tmp;
 
	for (dia = 0; dia < n; dia++) {
		max_row = dia, max = A(dia, dia);
 
		for (row = dia + 1; row < n; row++)
			if ((tmp = fabs(A(row, dia))) > max)
				max_row = row, max = tmp;
 
		swap_row(a, b, dia, max_row, n);
 
		for (row = dia + 1; row < n; row++) {
			tmp = A(row, dia) / A(dia, dia);
			for (col = dia+1; col < n; col++)
				A(row, col) -= tmp * A(dia, col);
			A(row, dia) = 0;
			b[row] -= tmp * b[dia];
		}
	}
	for (row = n - 1; row >= 0; row--) {
		tmp = b[row];
		for (j = n - 1; j > row; j--)
			tmp -= x[j] * A(row, j);
		x[row] = tmp / A(row, row);
	}
#undef A
}

int main(void)
{
	double a[] = {
		1, 1, -1, 
		2000, 0, 1000,
		0, 5000, 1000
	};
	double b[] = { 0, 5, 15 };
	double x[3];
	int i;
 
	gauss_eliminate(a, b, x, 3);
 
	for (i = 0; i < 3; i++)
		printf("%g\n", x[i]);
 
	return 0;
}

I run it:

0.000882353
0.00235294
0.00323529

See also: https://en.wikipedia.org/wiki/Gaussian_elimination, http://mathworld.wolfram.com/GaussianElimination.html.

But a fun with SMT solver is that we can solve these equations without any knowledge of linear algebra, matrices, Gaussian elimination and whatnot.

According to the source code of Z3, it can perform Gaussian Elimination, perhaps, whenever it can do so.


→ [list of blog posts]

Please drop me email about any bug(s) and suggestion(s): dennis(@)yurichev.com.