Transients vs Sympy

For variety, I switched from SageMath to Sympy since Sympy is a small more focused package. However, I’m struggling with it’s lack of seeming lack of equational reasoning. When you define an equation via Eq(lhs,rhs), it’s not a “first class” object – ie. you can’t divide it by 2 and expect both sides to be half as big, and you have to manually say “subs(titute) eq.lhs with eq.rhs” each time. Finally, I haven’t found good ways to say “substitute just this one occurrence of foo” – sympy wants to do them all. This is all totally anathema to how I think about maths – I’m used to pointing at some subexpression and saying “and use equation 2 to simplify that” etc.

To try out Sympy, I took up the challenge in Chapter 1 of Nahin’s Transients or Electrical Engineers book, in which he conjures up a second-order differential equation for a simple R/L/C ciruit and invites the reader to check he’s not lying. This was a decent workout, since 1) it involves first- and second-derivatives, 2) it involves plenty of equations and manipulating expressions. To be honest, I found this hard. I tried doing it on computer and got nowhere. I then did it by hand using paper and pen to figure out the steps. Then I went back to the computer and had to tackle several “how do I get the computer to do x” puzzles before I got anywhere.

In the end, I did indeed verify that Nahin is telling the truth (phew!) but it was hours and hour of headscratching work. It’s nice to have the computer verify I haven’t err’d in my derivation, but this was not “computer makes task easier” stuff – quite the opposite! I guess the nice bit comes later when I can get the computer to actually solve the differential equations for concrete cases. But this ‘abstract algebra manipulation’ stuff is painful.

In the end, the “proof” looks like..

from sympy import *
from import R,L,C,t

i = Function("i")(t)
i1 = Function("i1")(t)
i2 = Function("i2")(t)
v = Function("v")(t)
u = Function("u")(t)

# Chapter 1 of Transients For Electrical Engineers

# Equations which describe the given "R ser L par (R ser C)" circuit (ie. KVL/KCL)
eq1 = Eq(i, i1 + i2 )
eq2 = Eq(u, i*R + v)
eq3 = Eq(v, L * i1.diff())
eq4 = Eq(v, i2*R + (1/C) * Integral(i2,(t,0,t)))

# The book states a large differential equation with 'input' voltage u(t) + derivative on left,
# and current i(t) (+derivatives) on right
(u.diff(t,2) + (R/L) * u.diff() + 1/(L*C)*u)

# and challenges you to check the equality does actually follow from above 4 equations

# First step is easy; there's only one equation involving 'u' so let's apply that everywhere
_.subs( eq2.lhs, eq2.rhs).simplify()

# Now we have to be tactical; we have v and v' and v'' but we want to do different things to each
# I don't know how to tell sympy to "just substitutes the first instance, not all instances"
# So start by substituting for v'' (since that leaves v' and v alone for now) and use eq4 the capacitor eqn
_.subs(eq4.lhs.diff(t,2), eq4.rhs.diff(t,2)).simplify()

# And now we can tackle v' and then v - both use eq3 (the inductor one)
# This leaves us with just currents, albeit a mixture of i/i1/i2
_.subs(eq3.lhs.diff(t), eq3.rhs.diff(t)).simplify()

# Now we want to group + collapse the i1's and i2's using eq1.
# But I don't know how to get sympy to "see" that we have lots of (i1+i2) patterns 
# So as a hack, we can always subtract something then add something equal to it and preserve equality, right?
# So we use that to turn i1+i2 into i ...
(_ + diff(eq1.rewrite(Add),t)/C).expand()

# Getting there .. now we need to sweep up the twice-differentiated cases..
(_ + R*eq1.rewrite(Add).diff(t,2)).expand()
# w^5 - which was what we wanted

Sagemath vs Common Emitter Amplifier

Using Sagemath for electronics is nice because it forces you to be explicit about your assumptions and calculations. Ideally, you’d start by telling Sagemath about Maxwell’s laws and build up from there. For DC conditions, like calculating the bias for a common emitter amplifier we can assume no changing magnetic fields, so Faraday’s law of induction simplifies to Kirchhoff’s Voltage Law. We also need Ohm’s Law for the voltage-current relationship in resistors. I supplied the series/parallel resistance formula (but those too could come from KVL/KCL + Ohms law).

First let’s tell Sagemath about Kirchhoff’s Voltage Law:

from functools import reduce
from operator import add, neg

# Teach sagemath about physics (Faraday's law of induction in its KVL guise)

# Kirchhoff voltage law says sum of voltage drops around a loop is zero.
# Note: this is really just Faraday's law of inducion, for the special case of no changing magnetic fields
# Note: this is taking in symbolic expressions, and returning a symbolic expression
def KVL( voltages ): 
    return reduce(add,voltages) == 0

# Convenience case for where we have one voltage source, then we give remaining voltage drops as positive values 
# ie. KVL_vsource(5, [1,2,x] means a 5v source, then drop of 1v then drop of 2v then drop of x volts.
def KVL_vsource( vsource, vothers ):
     return KVL( [vsource] + list(map(neg,vothers)))

Next, let’s tell it about series/parallel resistors and Thevenin equivalent for voltage dividers (again, would be nice to derive this from KVL/KCL/Ohm)

# Kirchhoff's Current Law (aka conservation of charge) and Ohms Law give composition rules for impedance
# TODO: derive directly from KCL/Ohm.
recip(x) = 1 / x
ipar(a,b) = recip( recip(a) + recip(b))
iser(a,b) = a + b

# Thevenin equivalent voltage/resistance for divider (r1 and r2)
Vth_divider(v, r1, r2) = v * r2 / (r1+r2)
Rth_divider(r1, r2) = ipar(r1,r2)

Now we can pick some values for a common-emitter amplifier

# Transistor params
beta = 100
Vbe = 0.7

# Common emitter amplifier, biased using divider, with emitter feedback
Vcc = 5
Rtop,Rbot = 10e3, 3.3e3
Re = 470
Rc = 2e3

And calculate the quiescent voltages + currents:

# Turn the divider into thevenin equivalent
Vth = Vth_divider(Vcc, Rtop, Rbot )
Rth = Rth_divider(Rtop, Rbot)

# To get Ib, we use Kirchhoff's Voltage Law for loop around base/emitter, and solve for Ib
Ib = var("Ib")
eqn = KVL_vsource( Vth, [ Ib*Rth, Vbe, beta*Ib*Re])
print("KVL equation is ", eqn)
Ib = solve(eqn, Ib, solution_dict=True)[0][Ib].n()
print("So Ib is", Ib*1000, "mA")
# Finally, get voltage drop over Rc to get output voltage.
Ic = beta * Ib

Vc = Vcc - Ic * Rc
print("Vc is ", Vc)

vce_eqn = KVL_vsource(Vcc, [Ic*Rc, Vce, Ic*Re])
Vce = solve(vce_eqn, Vce, solution_dict=True)[0][Vce].n()
print("Vce is ", Vce)

Here we’re just using Sagemath as a numeric calculator. It’s kinda neat to package up Kirchhoff’s Voltage Law so we can just state all the knowns and let Sagemath find the unknowns, without us having to do algebra. But fundamentally, this is stuff we could easily do in excel/python/etc.

But with Sagemath we can stick with symbolic expression. Let’s start by keeping the potential divider as-is, but we’ll make Rc and Re symbols.

# Common emitter amplifier, biased using divider, with emitter feedback
Vcc = 5
Rtop = 10e3
Rbot = 3.3e3
Re = var("Re")
Rc = var("Rc")

# Turn the divider into thevenin equivalent
# Note: some textbooks just ignore supply impedance, ie. assume it's a perfect voltage source which has fixed voltage
# regardless of how much current is drawn; this isn't the case for resistor-divider.  We can either make heuristic assumptions
# (eg. assume load has "high enough" impedance, so draws little current) or we can calculate exactly via Thevenin equivalent.
Vth = Vth_divider(Vcc, Rtop, Rbot )
Rth = Rth_divider(Rtop, Rbot)

# To get Ib, we use Kirchhoff's Voltage Law for loop around base/emitter, and solve for Ib
Ib = var("Ib")
eqn = KVL_vsource( Vth, [ Ib*Rth, Vbe, beta*Ib*Re])
Ib = solve(eqn, Ib, solution_dict=True)[0][Ib]
Ic = beta * Ib

Vc = Vcc - Ic * Rc
print("Vc is ", Vc)
solve([Vc == Vcc/2, Ic*Re==Vcc/10], [Rc,Re], solution_dict=True)

# Output:
Vc is  -719/10*Rc/(133*Re + 3300) + 5

So now the collector voltage is left as an expression, dependent on the (as yet unchosen values of Rc and Re).

Rather than picking them manually, we can express constraints:

solve([Vc == Vcc/2, Ic*Re==Vcc/10], [Rc,Re], solution_dict=True

# Output
[{Rc: 13750/9, Re: 2750/9}]

So we’ve effectively just stated the design rules of “collector should be midrail” and “Ve should be about 10% of supply”, and Sagemath has found the values of Rc and Re which “make it so”. Nice.

How about if we wanted IC=20mA and Vc at midrail?

solve([Vc == Vcc/2,  Ic==0.01], [Rc,Re], solution_dict=True)

# Output
[{Rc: 250, Re: 3890/133}]