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 sympy.abc 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()
_.subs(eq3.lhs,eq3.rhs).expand().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
```