### Solving simultaneous equations by perturbation

Rupert Swarbrick <rswarbrick <at> gmail.com>

2014-08-28 21:03:08 GMT

This is just an email to say "Hey, that's easier than I
expected...". I've been messing around with some electronics and had
some ideas about making a current monitor to go with my (ancient) bench
power supply.
It turns out you really want several decades of resistors: in my case, I
wanted to measure low currents (100uA) up to reasonably high ones
(3A). I decided to use six shunt resistors, with values
0.004, 0.4, 4, 40, 400, 4000 ohms.
I'll ignore the problem of switching in different resistors, except to
say that you don't want to include the FET switches in the voltage
you're measuring because you don't have much control over the on
resistance of the transistor. So I came up with the topology that I've
sketched in the attached pdf (I only drew 3 resistors because I'm lazy,
and I threw it together quickly in Inkscape so please don't judge my
artistic skills too harshly!)
After staring at scribbled circuits on pieces of paper for a minute I
realised a problem: I don't want ridiculously large vertical resistors
(because of noise problems). Let's say the vertical resistors are each
maybe 100e3 ohms. But then suppose I'm using the 4k resistor: it's now
basically in parallel with five other 200k resistors, which will drop
the measured resistance significantly.
At this point, you might be wondering if I've sent this email to the
right list but, yes, that bit comes next!
I realised I could write down six simultaneous equations for r1
(nominally 0.004 ohm) to r6 (nominally 4000 ohm) and solve them. After a
couple of false starts, I wrote this code to generate the simultaneous
equations:
par(x,y) := x*y/(x+y);
partitions_1n (lst) :=
block([acc: [], before: []],
while not(emptyp(lst)) do (
acc: cons([first(lst),
append(reverse(before), rest(lst))], acc),
before: cons(first(lst), before),
lst: rest(lst)),
reverse(acc))$
res_net_eqn (pair, stopper) :=
block([nominal: first(pair), lh, rh],
lh:
if not (emptyp (second (pair))) then
ratsimp(par(1/lsum(1/(2*stopper + first(eqn)), eqn, second(pair)),
first(nominal)))
else
first(nominal),
rh: second(nominal),
num(lh) = rh * denom(lh))$
res_net_eqns (nominals, stopper) :=
makelist(res_net_eqn(pair, stopper),
pair, partitions_1n(nominals));
eqns: res_net_eqns([r1=4/100,r2=4/10,r3=4,r4=40,r5=400,r6=4000], 1/G)$
(For an explanation of partitions_1n, skip to the end)
Here "G" is the conductance of one vertical resistor. The resulting
equations are a bit unpleasant. For example:
(%i44) first(eqns);
(%o44) r1*r2*r3*r4*r5*r6*G^5+((((2*r1*r3+2*r1*r2)*r4+2*r1*r2*r3)*r5
+2*r1*r2*r3*r4)
*r6
+2*r1*r2*r3*r4*r5)
*G^4
+(((4*r1*r4+4*r1*r3+4*r1*r2)*r5
+(4*r1*r3+4*r1*r2)*r4+4*r1*r2*r3)
*r6
+((4*r1*r3+4*r1*r2)*r4+4*r1*r2*r3)*r5
+4*r1*r2*r3*r4)
*G^3
+((8*r1*r5+8*r1*r4+8*r1*r3+8*r1*r2)*r6
+(8*r1*r4+8*r1*r3+8*r1*r2)*r5
+(8*r1*r3+8*r1*r2)*r4+8*r1*r2*r3)
*G^2
+(16*r1*r6+16*r1*r5+16*r1*r4+16*r1*r3+16*r1*r2)*G
+32*r1
= ((((((r2+r1)*r3+r1*r2)*r4+r1*r2*r3)*r5+r1*r2*r3*r4)*r6+r1*r2*r3*r4*r5)*G^5
+((((2*r3+2*r2+4*r1)*r4+(2*r2+4*r1)*r3+4*r1*r2)*r5
+((2*r2+4*r1)*r3+4*r1*r2)*r4+4*r1*r2*r3)
*r6
+(((2*r2+4*r1)*r3+4*r1*r2)*r4+4*r1*r2*r3)*r5+4*r1*r2*r3*r4)
*G^4
+(((4*r4+4*r3+4*r2+12*r1)*r5+(4*r3+4*r2+12*r1)*r4+(4*r2+12*r1)*r3+12*r1*r2)
*r6
+((4*r3+4*r2+12*r1)*r4+(4*r2+12*r1)*r3+12*r1*r2)*r5
+((4*r2+12*r1)*r3+12*r1*r2)*r4+12*r1*r2*r3)
*G^3
+((8*r5+8*r4+8*r3+8*r2+32*r1)*r6+(8*r4+8*r3+8*r2+32*r1)*r5
+(8*r3+8*r2+32*r1)*r4+(8*r2+32*r1)*r3
+32*r1*r2)
*G^2+(16*r6+16*r5+16*r4+16*r3+16*r2+80*r1)*G+32)
/25
So I wasn't that surprised when shoving them straight into Maxima's
solve() function wasn't particularly successful - I managed to warm the
room up, but didn't achieve much else...
(If you're a proper analyst or physicist, you might want to avert your
eyes at this point: I'm far more comfortable doing algebra, so this
might be painful to watch)
Of course, then I realised that G is small. If I were a proper
physicist, I'd rescale to make it a dimensionless quantity, but I'm
feeling lazy and I've been at work all day. So, hmm. Perturbation
theory. Can't be that hard, right?
And with Maxima it turns out this problem is quite easy!
simple_soln_p (soln, vars) :=
is ((length (soln) = 1) and
setify(map(lhs, first(soln))) = setify(vars))$
flip_lists (lsts) := args (transpose (apply (matrix, lsts)))$
taylor_stratify_one (expr, x, x0, order) :=
block ([expansion: taylor(expr, x, x0, order), pieces],
pieces: makelist (x^n*coeff (expansion, x, n), n, 0, order),
if not (equal (expansion, apply("+", pieces)))
then
error ("Couldn't stratify expression")
else subst(1, x, pieces))$
taylor_stratify_eqns (eqns, x, x0, order) :=
flip_lists (makelist (taylor_stratify_one(eqn, x, x0, order), eqn, eqns))$
perturb_solve(eqns, vars, eps, order) :=
block([newterms: makelist([v,
makelist(gensym()*eps^k, k, 0, order)], v, vars),
newvars, newvar_subst, neweqs],
newvar_subst: makelist(first(pair) = apply("+", second(pair)),
pair, newterms),
neweqs: subst(newvar_subst, eqns),
/* newvars are grouped by order, not original variable */
newvars: flip_lists (subst(1, eps, map (second, newterms))),
block([strata: taylor_stratify_eqns (neweqs, eps, 0, order),
known: []],
for pair in flip_lists([strata, newvars])
do
block([stratum, these_vars, rewritten, soln],
[stratum, these_vars]: pair,
rewritten: subst(known, stratum),
soln: solve(rewritten, these_vars),
if not (simple_soln_p (soln, these_vars)) then
error("Can't solve ", rewritten, "for variables ", these_vars,
". Got: ", soln),
known: append(known, first(soln))),
subst(known, newvar_subst)))$
The "stratify" stuff uses taylor() to expand an expression in powers of
some variable, then collects coefficients for the different orders. Then
perturb_solve expands the variables in eps, stratifies the resulting
equation in terms of powers of eps and then solves a layer at a time,
From order zero upwards. I suspect there are equations where this
strategy fails miserably, but it turns out this isn't one of them.
And now to actually run the code:
(%i48) subst(1/100e3, G,
perturb_solve (eqns, [r1, r2, r3, r4, r5, r6], G, 1));
Evaluation took 0.0720 seconds (0.0730 elapsed) using 28.288 MB.
(%o48) [r1 = 0.04000004,r2 = 0.400004,r3 = 4.0004,r4 = 40.04,r5 = 404.0,r6
= 4400.0]
"First order?", I hear you scoff. So:
(%i49) subst(1/100e3, G,
perturb_solve (eqns, [r1, r2, r3, r4, r5, r6], G, 5));
Evaluation took 0.6480 seconds (0.6470 elapsed) using 295.398 MB.
(%o49) [r1 = 0.04000003980823639,r2 = 0.400003980860735,
r3 = 4.000398123173552,r4 = 40.03984945656259,r5 = 404.0224806732696,
r6 = 4444.219040588743]
Awesome!
Anyway, if anyone is still reading
(1) From bitter experience, I know that programming in Maxima is rather
difficult if you're trying to make robust code. But just hacking
something together is pretty easy. I wrote all the stuff above in
about 2 hours. (And was making up the maths as I went along:
someone who actually knew what they are doing would have been much
quicker)
(2) When writing code like this, my preferred approach is to sprinkle
print() statements liberally and then run it on small
examples. Most mistakes will give print-outs that are "clearly
wrong". (You can delete the print calls before you show anyone else
and pretend it worked first time...)
(3) The "flip_lists" function is kind of cute. I started writing a
zip() function and suddenly realised it's just matrix
transposition!
(4) It's worth indenting the code carefully. Maxima's syntax eschews
props for the weak-minded like curly braces (Bah! C++! Trivial!)
and has rather unhelpful syntax error messages when you get it
wrong. So make it as difficult to get wrong as possible.
(5) My partitions_1n function is really ugly. What it does is give you
all the ways of partititioning a list into "1 and n". That is
(%i50) partitions_1n([1,2,3,4]);
(%o50) [[1,[2,3,4]],[2,[1,3,4]],[3,[1,2,4]],[4,[1,2,3]]]
I couldn't think of a better way to express it than this
though. Anyone got a prettier version?
I hope this post was interesting / amusing for someone. The code was
definitely fun to write, and it answered a question I actually wanted to
answer. Woohoo!
Rupert

------------------------------------------------------------------------------
Slashdot TV.
Video for Nerds. Stuff that matters.
http://tv.slashdot.org/

_______________________________________________
Maxima-discuss mailing list
Maxima-discuss <at> lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/maxima-discuss