Symbolic Algebra in Python¶

Computers are inherently excellent at doing numeric calculations. By default, they are not nearly as capable of doing the types of symbolic manipulation you might do in an algebra or calculus class, and thus software like Mathematica and Maple exists. We can get similar functionality in Python (for free!) by using the Sympy module. This notebook is intended as a sort of crash course in the basics of Sympy functionality.

Setup ¶

If you installed Python through the Anaconda distribution, you should already have Sympy on your system and you just need to import it:

import sympy as sym

I almost always import sympy as sym, but if you are only going to be using Sympy in a notebook, you can eliminate some repetitive typing by using from sympy import * instead. Use this only if you understand name-spaces enough to keep yourself out of trouble though!

The above is all we'd need to get Sympy working. However, we can also ask it to use LaTeX in our notebook to output the equations all pretty! To do so, we just need to run

sym.init_printing()

Now we are good to go! Sympy isn't the lightest of modules, so I'd generally recommend you only import it if you are actually going to be using it in a notebook.

Jupyter notebooks have a built in help functionality through the use of ?. If you ever need to get more information about what inputs a function has or some examples on how to use it, a good place to start is always with function_name?. This will pop up a help window in the bottom of your notebook that you can use as a reference. For quick questions, it is almost always faster than googling for an answer!

Creating and Manipulating Functions ¶

When dealing with symbols in a programming environment, we immediately run upon a small conundrum. In programming we use variables to refer to various numeric values. For instance we might assign x = 10. However, when we are doing symbolic math, our "variables" are symbols, and we won't be assigning a value to them. In order to be clear about what type of variable we are creating, we have to specifically state what variables we want to be symbolic variables in Sympy. We do so by using sym.symbols:

x, y, z, t = sym.symbols('x, y, z, t')

The values on the left are what the programming variables will be named. The values in quotes on the right are the symbolic variables those programming variables represent.

Note that they don't necessarily have to match, so if you wanted to make your life miserable, you could do something like x, y, z = sym.symbols('z,x,y'), at which point whenever you used the variable x, you'd be referring to the symbolic variable z. I don't recommend this approach...

Creating Expressions¶

Once our variables are defined, we can use them to create functions!

3*x + (4*y)/(z**2)

You can see the pretty printing at work in the above expression. Note that you still need to use things like ** for "raised to the power of" operators. Also, if you are needing other special functions like sin or sqrt, you'll need to use the Sympy versions. Below we'll assign our function to the variable expr.

expr = 3*x + sym.cos(4*y)/sym.sqrt(5*z)
expr   # if you want the function to display after you create it, you need this

Combining and Manipulating Expressions¶

Once you've defined expressions, you can manipulate them just like other variables. So you can add, subtract, multiply, etc. For example:

combined = expr**2 + 2*expr
combined

There are however some particular features of equations that we might want to take advantage of. One of the first, and most useful, is the simplify function. Simplify tries its best to rewrite the expression you've given it into the most simple form. It does a decent job, though sometimes "simple" is subjective.

combined.simplify()

We can also expand an expression into all its parts, which will "foil" everything out:

combined.expand()

Or, in the opposite direction, we can factor an expression to try to factor out as many mutual terms as possible. Here note that I called sym.factor instead of using .factor(). You can generally use whichever notation you prefer! (However, sym.factor? will give you much more useful help than combined.factor?.)

sym.factor(combined)

Solving Equations¶

Perhaps one of the nicest aspects of Sympy is its ability to numerically solve an expression for a particular variable, or to solve systems of equations for several variables. The notation is very similar to Mathematica, where you provide first the expression and then the desired variable(s) to solve for.

newexpr = x**4 - 5*x**2 + 2
sym.solve(newexpr, x)

Two things to note:

Solutions are returned as a list, so you could save them as a variable and then access them with an index later
I never set what the expression was equal to! By default, Sympy assumes expressions are equal to 0.

If you have an expression that is not equal to 0, you can either rearrange it to get everything on one side or you can use sym.Eq( lhs, rhs ) where the two arguments are the left and right hand sides of your equation. For example, to get the 3rd solution to $x^4 - 5x^2 + 2 = -3$ I could go:

sols = sym.solve(sym.Eq(newexpr,-3), x)
sols[2] # 2 for the 3rd solution since Python starts at 0

Making Substitutions¶

Finally, a basic task that you often might like to do is substitute different values in for your symbolic variables. Maybe you want to evaluate your expression when $x=3$, or maybe you want to make a change of variables to $x=t^2$. For these situations, the subs command is Sympy is your friend. Suppose I have an expression $x^2 + 4x - 1$:

expr3 = x**2 + 4*x - y
expr3

I can substitute in a value like:

expr3.subs(x,3) # variable first, then value

Or you can substitute multiple values or terms like:

expr3.subs( [(x,t**2), (y,10)]) # sub x -> t**2 and y -> 10

# If you like, you can also use a dictionary for this and it is a bit cleaner imho
# expr3.subs({x:t**2, y:10})

Basic Calculus ¶

In addition to make algebra, frequently we might want a little help or desire to check our work on some calculus. Sympy is fully capable of doing derivatives and integration given the frameworks we have created. Derivatives utilize the diff command, and can again be written in either sym.diff() or expr.diff() formats. Recalling that our 3rd sample expression looked like:

expr3

Derivatives¶

We can take a derivative with respect to $x$:

expr3.diff(x)

Note that this is really the partial derivative, as it treats $y$ like a constant! If your function has only a single variable, you don't need to specify the derivative variable. You can do multiple differentiations easily by passing a number after the variable. So, for example, the 2nd derivative of the above would just be:

sym.diff(expr3, x, 2)
# or expr3.diff(x,2)

Integration¶

We can do integration in a similar and easy fashion. For an indefinite integral:

expr3.integrate(y)

If you want limits along with your integral, include the variable of integration and the upper and lower limits in a tuple. So $$ \int_0^{10} x^2 + 4x-y \, dy$$ is:

expr3.integrate( (y,0,10) )

We can even do multiple integration steps easily be chaining together more integration variables.

So $$ \int_0^{10} \int_{-2}^2 x^2 + 4x - y \, dx \, dy $$ will just be:

expr3.integrate( (x,-2,2), (y,0,10) )

Vector Calculus ¶

Finally, we are going to be doing some vector calculus this semester, and it would be nice to have a way to compute or check our work when finding gradients, divergences and curls! Sympy still has us covered here, but it is going to be a little more convoluted (which, I suppose, makes since given that vector calc is just another layer atop normal calculus...). To begin, we need to import a new sub-module of sympy: the vector module.

import sympy.vector as sv

Now, when dealing with vectors, our coordinate system is extremely important. Thus everything vector related in Sympy has to start out with defining a coordinate system. For a standard Cartesion coordinate system, we can do:

R = sv.CoordSys3D('R')

If you are using an older version of Sympy, you might not have the CoordSys3D function. In this case you'll need to use CoordSysCartesian (I think) to get the old version. Then everything will mostly be the same except for when taking gradients, divergences and curls you'll need to specify the coordinate system as well. Using ? will help you get the syntax correct. Or you could just update Sympy.

This creates our Cartesian coordinate system and names it R. It also specifies that we want our unit vectors to display as $\hat{x}$, $\hat{y}$, etc. If we want to refer to unit vectors in this coordinate system, we will use R.i, R.j, and R.k (unfortunately we can't change this). So I could define a simple vector as:

v1 = 1*R.i + 2*R.j + 3*R.k
v1

And we see it shows up with the nice notation. If we want to refer to the variable $x$ in the $R$ coordinate system, we need to use R.x. So the vector $$ \vec{v} = x \hat{x} + y\hat{y} + z\hat{z} $$ would be written as:

v2 = R.x * R.i + R.y * R.j + R.z * R.k
v2

It is easy to get mixed up with using R.x for the $\hat{x}$ unit vector, so be careful. Even if you have an expression with standard Sympy symbolic $x$ or $y$'s in it, you need to convert (or substitute) those to R.x's and R.y's before trying to do any vector operations on it!

Now that you can create vectors, you can add, subtract, dot or cross vectors like you'd expect:

v1+v2

v1.dot(v2)

v1.cross(v2)

v1.magnitude()

Complex Vector Calculus¶

We can also do more complex vector calculus involving gradients, divergences and curls. Say we have the potential energy equation $$ U(x,y) = 5x + \cos(xy) $$ And we'd like to find an expression for the force. Doing so requires us to take a gradient, since $$ \vec{F} = -\nabla U $$ We'll go ahead and define our potential energy first and make sure it looks ok:

U = 5*R.x + sym.cos(R.x * R.y)
U

And then we can just use the vector modules gradient command to compute the gradient:

F = -sv.gradient(U)
F.simplify()

Now that we have a vector expression we can go ahead and calculate the divergence and curl of $\vec{F}$:

sv.divergence(F)

sv.curl(F) # Conservative forces yo!

Calculating Line Integrals¶

In a variety of circumstances, we might want to use Sympy to calculate a line integral. More recent versions of Sympy finally added a convenience method to make this trivially simple for most situations, but when needed you can still fall back on the older methods, so I'll include them below as well. Let us take an example where we have a force field that looks like $\vec{F}=3y\hat{x} + (3x-2y)\hat{y}$.

F = 3*R.y*R.i + (3*R.x-2*R.y)*R.j

The Easy Way¶

To find the scalar potential relative to the origin, we can simply do:

U = -sv.scalar_potential(F, R)
U

Note that this will only work for conservative forces, because you clearly never specified any sort of path to integrate along. It also assumed you wanted the origin to be your 0 point. If those situations do not apply to you, you will have to use the more involved method.

The Hard(er) Way¶

Here we need to be much more explicit about what paths we are following. One potential path we could follow would be that of $$y=x$$ If we wanted to end up at some generic point $(x',y')$, e can parameterize this easily by imagining that we have to travel for the origin to $(x',y')$ in some time $t$, so that: \begin{align} x &= 0 + x' t\\ y &= 0 + y' t \end{align} We can thus define the line we want to integrate over by:

line = R.x*t*R.i + R.y*t*R.j

Where I am dropping the primes as I write them here. If we wanted to end up at the arbitrary point $(x,y)$ then, we would just let $t$ run from 0 to 1.

We then need to evaluate $F$ on this line, which basically means substituting the $x$ and $y$ parameterizations into our force equation:

FonLine = F.subs({R.x:line.dot(R.i), R.y:line.dot(R.j)})

where I am taking the dot products to grab first the x-component, and then the y-component and substitute them in properly. It can help to look at what FonLine looks like before moving on:

FonLine

Finally, we want to do the integration. To get $d\vec r$, we'll do a quick derivative:

dr = sym.diff(line,t)
dr

And now we do the integration:

-sym.integrate(FonLine.dot(dr), (t,0,1))

And we get the same expected value we got above!! Yay! This process can be wrapped in a function for easier use. Also note that, to do a piecewise parameterized function, you'll need to break it into pieces. For instance, to integrate first in the horizontal, and then in the vertical direction we might do something like:

horizontal = (R.x*t)*R.i # x is changing here, y is 0
dl_h = sym.diff(horizontal,t)
part_h = sym.integrate( F.subs({R.x:horizontal.dot(R.i), R.y:horizontal.dot(R.j)}).dot(dl_h), (t,0,1))

vertical = R.x*R.i + R.y *t*R.j # constant x now, but y is changing
dl_v = sym.diff(vertical,t)
part_v = sym.integrate( F.subs({ R.x:vertical.dot(R.i), R.y:vertical.dot(R.j)}).dot(dl_v), (t,0,1))

U = -(part_h+part_v)
U

And we get the same thing again (as we should, since this force is conservative). If you wanted to use a different point as your 0 point, you'd just have to change up the parametrization.

Taylor Series¶

We'll come back to look at certain Sympy functions more when the classwork demands them, but of particular interest to us at this moment is how we can find Taylor series with Sympy. One dimensional Taylor series are simple:

expr = 6*x*sym.cos(x)
expr.series(x)

By default Sympy displays the resulting series up to the 6th order and about the point 0. If you want to specify the point explicitly or the max order:

expr.series(x,0,n=4)

Often times you might want to use the expression without the clumsy O term. Getting rid of it is rather intuitive:

expr.series(x,0,n=4).removeO()

Now for the unfortunate part. Currently, Sympy does not recognize how to deal with multi-variable Taylor expansions. This is unfortunate, because multi-variable Taylor expansions are ugly to do by hand:

$$f(x,y) = f(a,b) + f_x(a,b)(x-a) + f_y(a,b)(y-b) + \frac{1}{2!}\left(f_{xx}(a,b)(x-a)^2 + 2f_{xy}(a,b)(x-a)(y-b) + f_{yy}(a,b)(y-b)^2\right) + \cdots.$$

We can still jury-rig Sympy to do the multivariable expansion for us, but it'll just be a smidge more involved and need you to do some checking at the end. Still easier than doing it by hand though!

Say we want to take a multivariable expansion about the point (0,1) of the expression:

$$ 6xy - \frac{x^2}{y}$$

We can get to the correct place by doing TWO taylor expansions, one in the $x$ and one in the $y$, and we'll keep all terms second order and less:

expr = 6*x*y - x**2/y
taylor = expr.series(x,0,n=3).removeO().series(y,1,n=3).removeO()
taylor

Note the need for the ```removeO()``` commands. If you don't get rid of the order symbol, the next series will choke on it!

This is pretty close to what we want. But we really only want to keep the mixed 3rd order terms at max. Because we did the two series independently, Sympy doesn't realize that, so we have to weed things out accordingly. Note that the first term is technically something along the lines of $x^2\cdot y^2$ which would be to the fourth order, and thus past our limit. As such, you need to manually remove it! This can sometimes be easier to see if you expand everything:

te = taylor.expand()
te

So we can see here that the first two terms have combined $x$ and $y$ powers greater than 2 and could be ignored. If you want to form an expression using only the desired terms, you can show the arguments of the expansion with:

te.args

And then decide that you want to keep the 0, and 3 terms. Careful, because it seems these are not ordered the same as displayed!

te_wanted = te.args[0]+te.args[3]
te_wanted

Fourier Series¶

We determined the needed integrals to take to find the coefficients of the Fourier series, but Sympy actually has a built in function that we can take advantage of to make the process even easier. Since we'll want to be plotting these, I'll go ahead and import matplotlib and numpy here (we already have Sympy imported and initialized)

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('seaborn-darkgrid')

Fourier Series of normal, continuous functions¶

I'll begin by defining a expression in Sympy. Recall that you need to define your variables in Sympy so that it knows what to keep symbolic! Initially I'll look at the function: $$ y(x) = x^2 $$

x = sym.symbols('x')
y = x**2

Then I'll find the fourier series of $y$ on some interval of $x$ (say from $x=-2$ to $x=2$.) This sometimes can take a bit, especially for more complex functions!

S = sym.fourier_series(y,(x,-2,2))

If we want to look at the first 3 terms of the series:

S.truncate(3)

Sympy Plotting¶

Interesting thing about Sympy. For really basic plotting, it actually has it's own plotter (though it relies on matplotlib for it's backend). Sympy's plotter is nice it that it works more like Maple or Mathematica (or WolframAlpha), but it is not as customizable as the full matplotlib. So it's great for quick and dirty plotting, not so great if you need a really nice pretty image. Here I could quickly look at the first few terms:

sym.plot(S.truncate(3), (x,-2,2));

Note that I can't do many nice things like add more plots to the graph, tweak the line settings, etc in the sympy plotter. But it works to give a quick and fast outlook, and even labels the axes (albeit in a ugly fashion...)!

Matplotlib Plotting of Sympy functions¶

Say instead that you wanted the full power of matplotlib, but you liked working with the sympy expressions. In this case you'll need to use lambdify to convert your sympy function to a more standard pythonic expression:

orig = sym.lambdify(x,y,'numpy')
fs = sym.lambdify(x,S.truncate(3), 'numpy')
xs = np.arange(-2,2,0.01)
plt.plot(xs,fs(xs), label='3 Term Fourier Series')
plt.plot(xs,orig(xs), label='Original')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

Lambdify takes three terms: the variable, the function, and then a modules option that tells it to convert its cosine/sine/sqrt/exp etc functions to numpy's versions.

Fourier Series of Piecewise Functions¶

Frequently we want to approximate piecewise functions as fourier series. As such, we need a method to define what a piecewise functions looks like! Take for example a sawtooth wave. It initially rises at some rate, then immediately drops to zero:

We can do something similar as follows:

p = sym.Piecewise( (x, x<0.5), (0,True) )
xs = np.arange(0,1,0.01)
plt.plot(xs, sym.lambdify(x,p,'numpy')(xs));

Note that I had to go ahead and lambdify p in order to be able to plot it. Sympy's plotter doesn't not seem to handle piecewise functions well. The piecewise function works by checking from left to right each point on the function and taking the first value that equals true. Thus, as long as $x<0.5$ in our example, then the function increases as $y(x)=x$. But as soon as $x>0.5$, then the function defaults to $y(x)=0$ (the True statement). You can add as many conditions as you need in your piecewise function, just remember to slowly work your way from left to right!

We can of course still take the fourier series of this function:

S2 = sym.fourier_series(p, (x,0,1))

Looking at the first few terms:

S2.truncate(3)

You can look at individual terms using the aptly named term:

S2.term(1)

And plotting the first few terms: (if you aren't using seaborn you'll need to set the colors yourself to classic values, or just delete the color parts and they'll all be the same!)

plt.figure(figsize=(12,8))

plt.subplot(221)
f = sym.lambdify(x,p,'numpy')
plt.plot(xs, f(xs))
plt.title('Original')

plt.subplot(222)
f = sym.lambdify(x, S2.truncate(3), 'numpy')
plt.plot(xs, f(xs), color='C1')
plt.title('Fourier Series: 3 terms')

plt.subplot(223)
f = sym.lambdify(x, S2.truncate(8), 'numpy')
plt.plot(xs, f(xs), color='C2')
plt.title('Fourier Series: 8 terms')

plt.subplot(224)
f = sym.lambdify(x, S2.truncate(30), 'numpy')
plt.plot(xs, f(xs), color='C3')
plt.title('Fourier Series: 30 terms')

plt.show()

Matrix Manipulations¶

As we enter into the final two chapters of the semester, we'll be seeing and using more and more matrices in our calculations. Often times, you can save yourself a load of algebra (and mistakes!) by doing some of these simple calculations in Python. Both Sympy and Numpy can do various matrix manipulations, but we'll focus mainly on Sympy.

Basic Definition and Operations¶

Defining a matrix is Sympy is rather simple:

A = sym.Matrix([[1,2,3],[4,5,6],[7,8,9]])
B = sym.Matrix([[2,4,6],[8,10,12],[14,16,18]])
A,B

We can add and subtract these just as you'd expect:

A-B

A+B

You can also do matrix multiplication in a simple fashion:

A*B

Or scale by different factors:

5*A

Advanced Operations¶

Other common operations you might want to do include traces, transposes, determinants or inverses.

The trace, or sum of the values on the diagonal:

A.trace()

The transpose flips the rows and columns of the matrix:

A.T

Determinants come up frequently:

sym.det(B)

You can't take the inverse of a matrix whose determinant is 0, so I'll define a new matrix which is invertible:

C = sym.Matrix([[4,1,0],[1,10,0],[0,0,5]])
sym.det(C)

C.inv()

Eigenpairs¶

Finding your eigenvalues and eigenvectors is deliciously easy:

C.eigenvals()

This format gives you first the eigenvalue and then the multiplicity, or how many eigenvalues were equal to that value. Here we got all unique values. If you want both the eigenvalues and the eigenvectors, then:

C.eigenvects()

For each pairing you'll get the eigenvalue, the multiplicity, and then the eigenvector corresponding to that eigenvalue. Note that eigenvectors can have variable scaling, so you might not always get the exact same values as you would by hand. However, the ratios between each value should be the same! You can always grab individual values via:

There *are* cases where solving the eigenvectors via numerics will give you very tiny imaginary contributions. Usually you can just set these to zero, but be aware!

evec2 = C.eigenvects()[1][2][0]
evec2

Which grabs the eigenvector corresponding to the second eigenvalue. You could then always normalize it by dividing by the magnitude (norm):

(evec2/evec2.norm())

That's ugly, so if you want it in numbers, remember evalf:

(evec2/evec2.norm()).evalf()

Dealing with Coupled Oscillation Modes¶

Say I take the instance where I have 3 blocks in a line with springs connecting them (2 springs total). I've worked out the matrix equations, and would like to visualize the different modes.

You can use either Sympy or Numpy to find your eigenvectors here, though I'll work mainly with Sympy. I'll begin by defining my K and M matrices:

K = sym.Matrix([[1,-1,0],[-1,2,-1],[0,-1,1]])
K

M = sym.Matrix([[1,0,0],[0,2,0],[0,0,1]])
M

If your M matrix is *not* equal to the identity matrix, then finding your normal frequencies and modes is *similar* but **not** identical to finding eigenvalues and vectors! So don't rely on programs to find the eigenvalues blindly!

Our $M$ matrix here is not equal to the identity, so we'll crunch the values ourselves. Let's form our required matrix:

ω = sym.symbols('ω', negative=False)
KmM = K - ω**2*M
KmM

As always, we know solutions will exist when the determinant is equal to 0. So finding the determinant and then solving the system for $\omega$:

sols = sym.solve(sym.det(KmM),ω)
sols

Note that by specifying that $\omega$ couldn't be negative above we avoided some negative values here. These are my values for $\omega$, so we have that $$\begin{align*} \omega_1 &= 0 \\ \omega_2 &= 1 \\ \omega_3 &= \sqrt{2} \end{align*} $$ Now we just need to solve for the correct corresponding $\vec{a}$ matrices! For $\omega_1$:

x,y,z = sym.symbols('x y z')
a = sym.solve(KmM.subs(ω,sols[0])*sym.Matrix([x,y,z]),(x,y,z))
a

Note that in the above I substituted in the first solution ($\omega=0$) for $\omega$ in the matrix equation, multiplied by a generic x,y,z vector, and then solved for when that would be equal to 0. We see that Sympy gives us a list of conditions when this will be true. I saved this list of conditions as a value for future use!

This lets us use the substitute command to set $x=z$ and $y=z$ in our matrix. Then I let $z=1$:

e1 = sym.Matrix([x,y,z]).subs(a).subs(z,1)
e1