Preface
In electronics engineering, from time to time you have to use standard formulas to characterize your circuits. To what extent you need to calculate all parameters most often depends on the requirement.
For example, consider the formula for the -3dB cutoff frequency of a 1st order RC lowpass filter:
f_c=\frac{1}{2\pi RC}Although this equation is fairly simple and most people won’t have any problem solving it for any particular variable in a few seconds, it can serve as a basic example on how to solve an equation symbolically.
One of the easiest ways of performing this task is to use SymPy, a Python library for symbolic mathematics.
An introductory example
First, we need to tell SymPy what the formula is:
from sympy import *
f, R, C = symbols("f_c R C")
# Remember not to name this 'f' (we already have a symbol named f!)
eq1 = Eq(f, 1 / (2 * pi * R * C))By using the SymPy variant of $\pi$ instead of numpy’s version, we can obtain higher numerical accuracy and the resulting equations will be formatted in a more human-readable way as they contain $\pi$ as symbol instead of
In case you get an error message like
ImportError: No module named 'sympy'
You can install the package using pip (recommended, but easy_install works too): pip install sympy
For details on the functions being used, refer to the excellent SymPy tutorial.
Now we can evaluate the original formula. As this is an equation and we want to compute the left-hand side ($\frac{1}{2πRC}$), we need to extract the right-hand side using .rhs and substitute the variables by actual values (for this example, we’ll use $1 k\Omega$ and $1 {\mu}F$). Last but not least, we use the evalf() function to evaluate the formula to a single floating point value
f_result = eq1.rhs.evalf(subs={C: 1e-6, R: 1e3})After running this codeblock, the value of f_result is about 159.155 — this is correct, as it can be easily verified using
The next step is to let SymPy solve the equation for R (in this example):
from sympy.solvers import solve
eq1_R = solve(eq1, R)[0] # [0]: Just take the first solutionNote that for more complex equations there might be more than one solution and you might or might not need to select the proper one. In this case there is only one solution: print(eq1_R) prints 1/(2*pi*C*f_c).
In frontends like , SymPy is also able to use various methods of pretty-printing. You can initialize the best method available for the current frontend by using
from sympy import init_printing
init_printing()In IPython, the Formula will be printed as LaTeX, with eq1_R resulting in $$f_c = \frac{1}{2πCf_c}$$
On the command line, an ASCII-based output will be generated:
1
─────────
2⋅π⋅C⋅f_cNow we can simply continue with substituting and evaluating the new expression. We will use the values from the previous block to evaluate the rearranged formula – it is expected that the result for R will be exactly 1000$\Omega$. solve(...)[0] already returns the right-hand side of the equation so we don’t need to use .rhs like above.
r_result = eq1_R.evalf(subs={C: 1e-6, f: f_result})
The value of r_result is exactly 1000.0\Omega. For more complex equations, numerical errors might yield even if SymPy tries to avoid them as best as possible. You can often bypass small numerical errors being returned by using the chop=true argument to evalf(). A detailed explanation of the mechanisms involved would exceed the scope of this article, therefore it is recommended to read the SymPy reference on numerical evaluation.
Lazifying SI units
Have you ever written a physics exam where you deduced the correct equation by matching the units involved to the expected output unit? Regardless, this is not neccessary any more.
SymPy can not only lazify your handling of symbolic equation but also integrated a unit system.
For example, you can calculate the unit of the equation listed above:
from sympy.physics.units import *
print(1 / (ohms * farads))This prints \frac{1}{s} which is correct considering the above equation yields a result in Hz.
SymPy can not only compute the correct deduced unit but is also able to deduce that is called Hz and
print(find_unit(1 / (ohms * farads))) # ['Hz', 'hz', 'hertz', 'frequency']We can also use SymPy’s knowledge of units to get a unit out of. We’ll use eq1 as listed in the previous section of this article.
Currently, SymPy does not support full calculation with units (see this issue). However, the current unit support can already be useful in order to sanity-check your equations. In the future I’ll try to integrate UliEngineering’s unit system with sympy.