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.<\/p>\n\n\n\n
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.<\/p>\n\n\n\n
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.<\/p>\n\n\n\n
In the end, the “proof” looks like..<\/p>\n\n\n\n
from sympy import *\nfrom sympy.abc import R,L,C,t\n\ni = Function(\"i\")(t)\ni1 = Function(\"i1\")(t)\ni2 = Function(\"i2\")(t)\nv = Function(\"v\")(t)\nu = Function(\"u\")(t)\n\n# Chapter 1 of Transients For Electrical Engineers\n\n# Equations which describe the given \"R ser L par (R ser C)\" circuit (ie. KVL\/KCL)\neq1 = Eq(i, i1 + i2 )\neq2 = Eq(u, i*R + v)\neq3 = Eq(v, L * i1.diff())\neq4 = Eq(v, i2*R + (1\/C) * Integral(i2,(t,0,t)))\n\n# The book states a large differential equation with 'input' voltage u(t) + derivative on left,\n# and current i(t) (+derivatives) on right\n(u.diff(t,2) + (R\/L) * u.diff() + 1\/(L*C)*u)\n\n# and challenges you to check the equality does actually follow from above 4 equations\n\n# First step is easy; there's only one equation involving 'u' so let's apply that everywhere\n_.subs( eq2.lhs, eq2.rhs).simplify()\n\n# Now we have to be tactical; we have v and v' and v'' but we want to do different things to each\n# I don't know how to tell sympy to \"just substitutes the first instance, not all instances\"\n# So start by substituting for v'' (since that leaves v' and v alone for now) and use eq4 the capacitor eqn\n_.subs(eq4.lhs.diff(t,2), eq4.rhs.diff(t,2)).simplify()\n\n# And now we can tackle v' and then v - both use eq3 (the inductor one)\n# This leaves us with just currents, albeit a mixture of i\/i1\/i2\n_.subs(eq3.lhs.diff(t), eq3.rhs.diff(t)).simplify()\n_.subs(eq3.lhs,eq3.rhs).expand().simplify()\n\n# Now we want to group + collapse the i1's and i2's using eq1.\n# But I don't know how to get sympy to \"see\" that we have lots of (i1+i2) patterns \n# So as a hack, we can always subtract something then add something equal to it and preserve equality, right?\n# So we use that to turn i1+i2 into i ...\n(_ + diff(eq1.rewrite(Add),t)\/C).expand()\n\n# Getting there .. now we need to sweep up the twice-differentiated cases..\n(_ + R*eq1.rewrite(Add).diff(t,2)).expand()\n# w^5 - which was what we wanted<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-1532","post","type-post","status-publish","format-standard","hentry","category-general"],"_links":{"self":[{"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/posts\/1532","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/comments?post=1532"}],"version-history":[{"count":4,"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/posts\/1532\/revisions"}],"predecessor-version":[{"id":1539,"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/posts\/1532\/revisions\/1539"}],"wp:attachment":[{"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/media?parent=1532"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/categories?post=1532"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.nobugs.org\/blog\/wp-json\/wp\/v2\/tags?post=1532"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}