View source for SymPy/Simultaneous Equations

==Introduction==
This page focuses on using SymPy to find both the symbolic and the numeric solutions to equations obtained from electric circuits.  

== Initialization ==
To use SymPy, you will need to import the module.  The module comes with the Anaconda distribution of Python, but if you have a different installation, you may need to get the module first.  Also, you will likely want to use the "pretty printing" capabilities of SymPy, which requires initializing printing.  These two steps can be done with:
<syntaxhighlight lang=python>
import sympy as sym

sym.init_printing()
</syntaxhighlight>

== Declaring Symbols ==
While there are ways to bring in certain "standard" symbols in SymPy, it may be clearest to define them yourself.  You need to declare something as a symbol before using it on the right side of an assignment; as usual, if you create a variable as a combination of symbols, Python will duck type it as a symbol too.  The <code>sym.var()</code> function accepts a space or comma-separated string of variables and will create symbolic entities based on the code in the string.  
For instance:
<syntaxhighlight lang=python>
sym.var('x, y')
</syntaxhighlight>
will create variables <code>x</code> and <code>y</code>

== Defining Equations ==
When it comes to solving things later, there are two different ways to define an equation:
* If the equation is equal to zero, you do not explicitly have to include the "=0" part; you can simply define a variable to hold on to the expression.  If $$x+y=0$$ is one of the equations you are looking to solve, you can use:
<syntaxhighlight lang=python>
eqn1 = x + y
</syntaxhighlight>

If you have an expression that evaluates to other than zero, you can use sym.Eq(left, right) for that.  If the other equations for which you are trying to solve is $$x-y=3$$, you can write that in code as:
<syntaxhighlight lang=python>
eqn2 = sym.Eq(x - y, 3)
</syntaxhighlight>
Of course, you could also rearrange the expression by moving the right side to the left, in which case:
<syntaxhighlight lang=python>
eqn2 = x - y - 3
</syntaxhighlight>
would work!


==Solving Equations With SymPy==
To solve a system of linear algebra equations, you can use the <code>sym.solve(eqns, vars)</code> command.  The system of equations should be in a list, set, or tuple as should the list of unknowns.  Note that even if you don't specifically "care" about the value of one or more of the unknowns, you must list them; if you want to solve the two equations above to find x, you still need to explicitly list y as a unknown or else SymPy will assume that it is a symbolic, but known, value, thus constraining the system.  Here is an example solving the two equations above:
<syntaxhighlight lang=python>
soln = sym.solve((eqn1, eqn2), (x, y))
</syntaxhighlight>
If you run this in the console, <code>soln</code> will print out pretty.  If you put all this in a script, <code>soln</code> will not show up at all...  You can put <code>print(soln)</code> in a script, but it will not be pretty.  To make it pretty requires a little work since the solution set is a dictionary with the variables as keys.  To help, here is a Python function you can add to your code:
<syntaxhighlight lang=python>
def printout(solution, variables, fp=False):
    for var in variables:
        sym.pretty_print(var)
        if fp:
            sym.pretty_print(solution[var].evalf(fp))
        else:
            sym.pretty_print(solution[var])
    print()
</syntaxhighlight>
This will go through the solutions and print the variable followed by the solution.  If you have substituted in numbers, set fp to True to do floating point evaluation on the answers.  If you define this function in the console and then type
<syntaxhighlight lang=python>
printout(soln, (x, y))
</syntaxhighlight>
you will see the variables and their values.  

When the solutions are numbers, this does not seem that useful, but imagine that you want to solve:
<center>$$
\begin{align*}
ax+by&=c\\
dx+ey&=f
\end{align*}$$</center>
symbolically?  You will need to add symbols for a through f first, then redefine the equations to include those symbols; here is a sample code:
<syntaxhighlight lang=python>
import sympy as sym

sym.init_printing()

def printout(solution, variables, fp=False):
    for var in variables:
        sym.pretty_print(var)
        if fp:
            sym.pretty_print(solution[var].evalf(fp))
        else:
            sym.pretty_print(solution[var])
    print()

x, y = sym.var('x, y')
a, b, c, d, e, f = sym.var('a, b, c, d, e, f')

eqn1 = a*x + b*y - c # shows moving right to left side
eqn2 = sym.Eq(d*x+e*y, f) # Shows using sym.Eq to give left and right

soln = sym.solve((eqn1, eqn2), (x, y))
print(soln)
printout(soln, (x, y))
</syntaxhighlight>
The printed version gives:
<syntaxhighlight lang=text>
{x: (-b*f + c*e)/(a*e - b*d), y: (a*f - c*d)/(a*e - b*d)}
</syntaxhighlight>
which is a bit hard to interpret; the printout version shows the two fractions as mathematical expressions.

== Substituting Values ==
SymPy has a subs method that works on a variable.  For instance, if you have:
<syntaxhighlight lang=python>
import sympy as sym

sym.init_printing()

sym.var('a, b, c, d, e, f')
a = b + c * d + e ** f
</syntaxhighlight>
You have created $$a=b+cd+e^f$$.  If you want to know what this would be if $$b=2$$, you could write:
<syntaxhighlight lang=python>
a.subs(b, 2)
</syntaxhighlight>
to see that - note that this does not change $$a$$.  If there are multiple substitutions, you can put them in an interable (i.e. list or tuple) of variable-value pairs; the following two lines will produce the same results:
<syntaxhighlight lang=python>
a.subs( ( (b,2),(c,3) ) )
a.subs( [ [b,2],[c,3] ] )
</syntaxhighlight>
If you have an iterable of variables to substitute for and an iterable of values, you can use the zip command to create those pairs for you:
<syntaxhighlight lang=python>
varstosub = [b, c]
valstosub = [2, 3]
zipsub = zip(varstosub, valstosub)
a.subs(zipsub)
</syntaxhighlight>
There will be times when you want to make the same substitions to several symbolic representations at once.
Given that, here is another handy function for making substitutions that takes a list of variables and a list of values:
<syntaxhighlight lang=python>
def makesubs(solution, vartosub, valtosub):
    subsol = solution.copy()
    sublist = list(zip(vartosub, valtosub))
    for var in solution:
        subsol[var] = subsol[var].subs(sublist)
    return subsol
</syntaxhighlight>

== Notes ==
For the printout and makesubs commands, if there is only one variable involved, the "list" of variables will be a single item and not an iterable, you may need to put [ ] around the variable "list" in the function call to get it to work.

<!--
== Examples ==
===Solving Circuit Variables===
[[File:SampleStarter.png|thumb]]
Each of the samples in this section will involve solving for the power
absorbed by resistor <math>R_{4}</math> in the circuit at right
by using a different method.  The circuit will be presented several times
along the way to demonstrate how to label it properly for the given
method.
<br clear=all>

====Node Voltage Method====
=====Choose Ground Node =====
[[file:SampleNVM1.png|thumbs|Circuit now shown with ground node.]]
First, choose a ground node.  In this case, the large bottom node is
as good as any:
=====Label Each Node=====

Next, circle or otherwise mark each node:
\begin{center}
\epsfig{file=./XMAPLE/SampleNVMnodes.ps, scale=0.5}
\end{center}
At this point, you can see that there are five total nodes.  One of
them is the ground node.  From ground, you can determine the node
voltage of the node at the positive terminal of $\E{v}{s}$ since the
voltage source determines the different in voltage between ground and
that node.  Otherwise, you will need to assign labels to the remaining
three nodes.  In this case, use $\E{v}{x}$, $\E{v}{y}$, and $\E{v}{z}$
as node voltages:
\begin{center}
\epsfig{file=./XMAPLE/SampleNVMnodeslab.ps, scale=0.8}
\end{center}
  Note that other unknowns could have been selected so long as they
  gave equations for each of the nodes.
=====Determine Equations=====
Finally, you can write the KCL equations.  In this case, there are
three unknowns so you need to pick three of the nodes.  To help
choose, remember that voltage sources are bad for KCL equations
because {\it any} amount of current can go through a voltage source.
Conveniently enough, that eliminates two of the five nodes from
consideration ($n_1$ and the ground node) leaving three available.
The KCL equations, written by determining the amount of current {\it
  leaving} the node through each branch that passes through the node
boundary, are:

Assume you have a circuit for which the KCL equations are:
<center><math>
\begin{align}
\mbox{KCL, n}_2&: &
\frac{v_{x}-v_{y}}{R_{1}}+
\frac{v_{x}-v_{s}}{R_{2}}+
\frac{v_{x}-v_{z}}{R_{3}}+
\frac{v_{x}-0}{R_{4}}&=0\\
\mbox{KCL, n}_3&: &
\frac{v_{y}-v_{x}}{R_{1}}-
i_{t}&=0\\
\mbox{KCL, n}_4&: &
\frac{v_{z}-v_{x}}{R_{3}}+
\frac{v_{z}-0}{R_{5}}+
i_{t}&=0
\end{align}
</math></center>
These three equations can be put into Maple ''as is'' - again, no
need to set them up as a matrix if you are using Maple.  Maple can
solve for the three unknowns, and the power absorbed by <math>R_{4}</math>
will be the voltage across the resistor (<math>R_{x}</math>)
squared divided by the resistance.


====Branch Current Method====
For the BCM, start be determining the number of branches.  In this
particular case, there are five branches.  One of them - the very top
one - has an independent current source so its current is known.  The
other four are yet to be determined:
\begin{center}
\epsfig{file=./XMAPLE/SampleBCMbranches.ps, scale=0.5}
\end{center}
Without more information, no more currents can be labeled so you will
need to create an unknown that will help determine other branch
currents.  Labeling the current through $\E{R}{5}$ as $\E{i}{y}$, for
example, will not only allow you to note what the current through
$\E{R}{5}$ is but also allow you to write an equation for the current
through $\E{R}{3}$ using KCL at the node common to  $\E{R}{3}$,
$\E{R}{5}$, and  $\E{i}{t}$:
\begin{center}
\epsfig{file=./XMAPLE/SampleBCMbranchesy.ps, scale=0.5}
\end{center}
This still leave two unknown currents, and there is as yet no way to
solve for either of them.  Since we eventually want to know the power
absorbed by $\E{R}{4}$, it makes sense to label its current so call
that $\E{i}{x}$.  Once that is done, you can use KCL at the ground
node to determine the current through the branch with voltage source
$\E{v}{s}$:
\begin{center}
\epsfig{file=./XMAPLE/SampleBCMbranchesxy.ps, scale=0.5}
\end{center}
You may also want to do a quick check at the node connecting all four
resistors to make sure KCL holds there.  Calculating the current
leaving that node yields:
\begin{align*}
-(\E{i}{t})-(\E{i}{x}+\E{i}{y})+(\E{i}{t}+\E{i}{y})+(\E{i}{x})&=0
\end{align*}

At this point, you have two unknowns and six loops to choose from. For
the current methods, you must avoid current sources which knocks out
three loops (the top right primary loop, the right superloop, and the
overall superdooperloop\footnote{Well...what would {\it you} call
  it?}).  That leaves three loops and thus two loop equations.  In
this case, each of the primary loops will have four elements while the
bottom superloop has four, so you may choose the primary loops:
\begin{center}
\epsfig{file=./XMAPLE/SampleBCMloops.ps, scale=0.8}
\end{center}
\begin{align*}
\mbox{KVL, l}_1&: &
-\E{v}{s}+
\E{R}{2}\left(\E{i}{x}+\E{i}{y}\right)+
\E{R}{4}\left(\E{i}{x}\right)&=0\\
\mbox{KVL, l}_2&: &
\E{R}{4}\left(-\E{i}{x}\right)+
\E{R}{3}\left(\E{i}{t}+\E{i}{y}\right)+
\E{R}{5}\left(\E{i}{y}\right)&=0
\end{align*}
These two equations can be put into Maple and Maple can
solve for the two unknowns.  In this case, the power absorbed by $\E{R}{4}$
will be the current through the resistor ($\E{i}{x}$)
squared multiplied by the resistance.
\pagebreak

====Mesh Current Method====
For the MCM, start by labeling the mesh currents in the primary loops.  
In this particular case, there are three primary loops.  As to avoid
confusion with the unknowns in the previous examples, they are labeled
with subscripts $a$, $b$, and $c$:
\begin{center}
\epsfig{file=./XMAPLE/SampleMCMcurrents.ps, scale=0.8}
\end{center}
Meshes $b$ and $c$ are ready for KVL.  Mesh $a$ has a
problem in that there is a current source.  This can be solved by
using an auxiliary equation based on the current source - looking at
that mesh, it should be clear that:
\begin{align*}
\E{i}{a}&=-\E{i}{t}
\end{align*}
so you may use that as one of your equations.  If a current source
happens to be between loops, you will end up using two auxiliary
equations - one will be the equation to define the current source in
terms of two mesh currents, and the other will be a superloop equation
that combines the two meshes touched by the current source.

At this point, you have two meshes left for which to write KVL equations:
\begin{align*}
\mbox{KVL, l}_b&: &
-\E{v}{s}+
\E{R}{2}\left(\E{i}{b}\right)+
\E{R}{4}\left(\E{i}{b}-\E{i}{c}\right)&=0\\
\mbox{KCL, l}_c&: &
\E{R}{4}\left(-\E{i}{b}+\E{i}{c}\right)+
\E{R}{3}\left(-\E{i}{a}+\E{i}{c}\right)+
\E{R}{5}\left(\E{i}{c}\right)&=0
\end{align*}
The three equations - the auxiliary equation plus the two KVLs - can
be put into Maple and Maple can 
solve for the three unknowns.  In this case, the power absorbed by $\E{R}{4}$
will be the current through the resistor (now, $\E{i}{b}-\E{i}{c}$)
squared multiplied by the resistance.

<!--
\section{Assignment}
Create three worksheets - one each to solve the equations generated by
the circuit above with the three methods.  At the end of each worksheet,
there should be a variable named {\tt pabsR4XXX} where {\tt XXX} is
the method used to generate the equations.  Use the {\tt simplify} and {\tt
  expand} commands on this final symbolic representation for the power
absorbed and make sure each of the three methods produce the same
answer.  Then use Maple's {\tt subs} command to determine the power
absorbed assuming $\E{v}{s}$=5 V, $\E{i}{t}$=3 mA, 
$\E{R}{1}$=1.2 k$\Omega$, 
$\E{R}{2}$=2.2 k$\Omega$, 
$\E{R}{3}$=3.6 k$\Omega$, 
$\E{R}{4}$=4.7 k$\Omega$, and
$\E{R}{5}$=5.6 k$\Omega$.  Display your final answer with 4
significant figures (i.e. use the {\tt evalf[4]} command).


% LocalWords:  NVM BCM MCM KCL KVL knowns supernode supernodes Vals MySoln OIT
% LocalWords:  evalf superloop Cramer's Maplesoft xmaple startup eqn
-->