In this article, we will see how to make graphs with SymPy and Matplotlib.
Usually when we talk about SymPy we focus on the part related to solving mathematical calculations. But many forget that the library has some basic resources for building graphs, which we will see in this article.
Creating an equation #
To illustrate, I will use the ideal gas law in this article. The famous PV=nRT seen in high school. Let’s check how to make graphs that show the dependence of the volume of an ideal gas on pressure and temperature.
import sympy
from sympy.plotting import plot, plot3d
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
# settings for better outputs in the article, can be ignored
sympy.init_printing(
use_latex="mathjax",
scale=1.0,
order="grlex",
forecolor="Black",
backcolor="White",
)
Just for reference, the version of SymPy in this article is the most recent at the time I write.
sympy.__version__
'1.12'
Remember that the gas constant R has several possible values depending on the unit system to be used. It is very likely that you have used the value “0.082” several times in your life because it is the value associated with the atmosphere (for pressure) and liters (for volume), which are common in most countries for these quantities. The temperature is expressed in Kelvin. In reality, this is an approximate value to facilitate the accounts made manually or in simple calculators. Here, we don’t need to approximate, so let’s put all the significant ones:
gas_constant = 0.082057366080960 # atm L / (mol K)
Actually, we should not approximate. Never make approximations in intermediate steps, this will only lead to errors propagated in your calculations. But that’s a subject for another article. Just take advantage of the fact that you have a computer that handles numbers very well and don’t make approximations. If you need a “well-behaved” number for an answer, round only at the end.
Returning focus, let’s now create symbols for each quantity and create our equation. If you need any clarification about this step, see the article on equations with SymPy:
P, V, T, R, n = sympy.symbols('P V T R n', posivite=True)
ideal_gas_law = sympy.Eq(P*V, n*R*T)
ideal_gas_law
\(\displaystyle P V = R T n\)
Making clear the role of each symbol:
- P: pressure in atmospheres (atm)
- V: volume in liters
- n: amount of substance in moles
- R: gas constant in atm.l/(mol.K)
- T: temperature in Kelvin
Let’s solve this equation for the volume V:
sympy.solve(ideal_gas_law, V)
\(\displaystyle \left[ \frac{R T n}{P}\right]\)
See that the result comes out as a list with a single element, which is the resulting expression of isolating V from the equation on the left side of the equality. Let’s associate such expression with a variable:
volume = sympy.solve(ideal_gas_law, V)[0] # index 0 to get the first element
volume
\(\displaystyle \frac{R T n}{P}\)
There are three variables in the expression. In other words, V=f(P, n,T), volume is a function of pressure, amount of substance and temperature. Let’s start simple, obtaining an expression for temperature and amount of substance constant. To do this, let’s substitute 1 mol and 25 °C (298.15 K) into the expression:
molar_volume_25_celsius = volume.subs({R: gas_constant, n: 1, T: 298.15})
molar_volume_25_celsius
\(\displaystyle \frac{24.4654036970382}{P}\)
Now we are in a position to make a graph showing how the volume V changes with variations in pressure P kept constant n and T at the values of 1 mol and 298.15 K, respectively.
Plots with SymPy #
At the beginning of the article, we saw that we imported the plot method from
SymPy. Let’s pass the expression obtained for the volume as an argument and
indicate which variable, in this case P:
plot(molar_volume_25_celsius, P)
<sympy.plotting.plot.Plot at 0x7f7cf6290800>
There are some interesting observations about the generated graph:
- Those who are more used to Python should have recognized the “face” of Matplotlib in the graph. And, indeed, SymPy uses Matplotlib as a backend for graphs
- Note how the axes were named. The abscissa (the “x axis”) was named P, the
variable indicated when calling the
plotmethod. The ordinate (“y axis”) was named f(P), indicating that this axis is a function of P. As it is a mathematical library, SymPy assumes that y = f(x) and makes this explicit in the graph. It is possible to change these labels, but it is important to note the default representation. - Another characteristic of a mathematical library, the graph is centered at the origin (0, 0). This is not a Matplotlib standard, it is a modification made by SymPy.
- But perhaps the most striking feature is the fact that there is a part of the
graph in negative values. Obviously, this does not make sense physically, since
it is about pressure and volume. In fact, when creating the symbols at the
beginning, we put
positive = True. So why did the graph appear in the 3rd quadrant?
We see in the documentation of the plot
method
that, by default, the domain value range considered is (-10, 10), as can be
seen in the graph. Therefore, regardless of what we ask SymPy to assume when
creating the symbols, the graph will be generated for the range requested by
Matplotlib. But this is a practical explanation.
In reality, we have to keep in mind that all these assumptions that we ask SymPy to make, such as that a certain symbol has only positive values, are only valid within the context of SymPy operations. SymPy is a symbolic abstraction layer over pure Python. When it interacts with other libraries, such as Matplotlib, there is no way to guarantee coherent behaviors with the assumptions. The assumptions are valid within SymPy operations and are important, as we saw in the article on logarithms. However, be careful with your expectations when interacting with other libraries.
Considering these observations, let’s restrict the domain to a range that makes more sense. To achieve this, just pass our variable and the minimum and maximum values in a tuple:
plot(molar_volume_25_celsius, (P, 0.1, 4))
<sympy.plotting.plot.Plot at 0x7f7cf629ef60>
We already have the characteristic profile of an ideal gas isotherm. If we still want to restrict the graph further, to a range of smaller volume values, just pass the desired minimum and maximum values for the y axis to the ylim argument in the form of a tuple:
plot(molar_volume_25_celsius, (P, 0.1, 4), ylim=(0, 80))
<sympy.plotting.plot.Plot at 0x7f7cf629d6a0>
And what if we want another curve on the graph? For example, for 100 °C (373.15
K)? Simple, just create the expression and pass it to the plot method:
molar_volume_100_celsius = volume.subs({R: gas_constant, n: 1, T: 373.15})
plot(molar_volume_25_celsius,
molar_volume_100_celsius,
(P, 0.1, 4),
ylim=(0, 80))
<sympy.plotting.plot.Plot at 0x7f7cf5b93ce0>
We see that a second curve was created, referring to the new temperature. We could continue adding curves by doing the procedure above and, of course, improve the presentation of the graph by identifying the curves. However, it will soon become evident that interacting with Matplotlib through SymPy is very bureaucratic and limited. For example, it doesn’t make much sense to create an expression for each curve if the intention is, for example, to present curves for 5 temperatures on the graph. And this is not a criticism of SymPy, it simply is not the focus of the library to deal with graphs, they have already done a lot in putting a simplified way of generating graphs. The point is that more advanced customizations become much simpler if we extract the data from the SymPy objects and pass them directly to Matplotlib.
Plots with Matplotlib #
Let’s unleash the unlimited power of Matplolib.
Thanks, Palpatine, let’s go to the Matplotlib side of the force… Let’s see how this works.
The first step is to request the plot method not to present the graph itself:
plot(molar_volume_25_celsius, (P, 0.1, 4), ylim=(0, 80), show=False)
<sympy.plotting.plot.Plot at 0x7f7cf5ba6e40>
Note that only a Plot object is generated at a memory address. Therefore, we
can associate such object with a variable. Let’s create a variable with the
super creative name of p:
p = plot(molar_volume_25_celsius, (P, 0.1, 4), ylim=(0, 80), show=False)
p
<sympy.plotting.plot.Plot at 0x7f7cf5c21d90>
We can order the graph to be displayed from this variable:
p.show()
Now, to comprehend the next steps, let’s understand in a very simplified way the composition of a Matplotlib graph. When you look at a graph like the one above, think of it as:
- a white rectangle, which is the base of the figure
- an axis system above this rectangle
- the axis system has the line element, the labels of each axis, the values presented on the axis and everything else we see
- if wanted, new axis systems can be added to the same base figure
- and, if wanted, new lines (or points) can be added to a given axis
This is basically how graphs are created with Matplotib. More details can be seen in the Matplotlib API documentation and with a more rigorous language. The intention here is to understand that we can decompose a graph into more fundamental elements (objects in programming).
Thus, the object associated with our variable p is, in fact, a container! It
has all these requirements to create a Matplotlib graph. Let’s see all the
public attributes and methods associated with this object:
[x for x in dir(p) if not x.startswith('_')]
['annotations',
'append',
'aspect_ratio',
'autoscale',
'axis',
'axis_center',
'backend',
'extend',
'fill',
'legend',
'margin',
'markers',
'rectangles',
'save',
'show',
'size',
'title',
'xlabel',
'xlim',
'xscale',
'ylabel',
'ylim',
'yscale',
'zlabel']
Now, it’s all there. Information about the axes, labels, (almost) everything we talked about earlier. For example, let’s check the label of the y axis and the limit we requested (from 0 to 80):
p.ylabel, p.ylim
\(\displaystyle \left( f{\left(P \right)}, \ \left( 0.0, \ 80.0\right)\right)\)
Great! Now, there are two “suspicious” methods in this list: append and extend. They are two characteristic methods of lists, associated with the ability to increase a list by adding more elements. So… the object p was created with a list structure. And what will be the content of this list? Let’s see the first position of the list:
p[0]
<sympy.plotting.plot.LineOver1DRangeSeries at 0x7f7cf5c235c0>
Hooray! A LineOver1DRangeSeries object, in other words, with all the letters,
a line over a range of one-dimensional values. Translating, our curve. And, as
there is only one curve, it doesn’t make sense to have an element in the next
position, right?
p[1]
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
Cell In[25], line 1
----> 1 p[1]
File ~/.virtualenvs/chemistry/lib/python3.12/site-packages/sympy/plotting/plot.py:265, in Plot.__getitem__(self, index)
264 def __getitem__(self, index):
--> 265 return self._series[index]
IndexError: list index out of range
Right! So we already know that new curves are new elements in this list that is
the object p. And what are the public attributes and methods available for
each element of this list?
[x for x in dir(p[0]) if not x.startswith('_')]
['adaptive',
'depth',
'end',
'expr',
'get_color_array',
'get_data',
'get_points',
'get_segments',
'is_2Dline',
'is_3D',
'is_3Dline',
'is_3Dsurface',
'is_contour',
'is_implicit',
'is_line',
'is_parametric',
'label',
'line_color',
'nb_of_points',
'only_integers',
'start',
'steps',
'var',
'xscale']
We see, by the names, that they are information about the curve itself: its
color, number of points and the points themselves. And it is the points that we
want, to be able to build our own graphs directly with Matplotlib. We can get
the points with the get_points method (suggestive, isn’t it?):
# the number of values is quite long, so I used a little bit of NumPy here
# just to show the first three and the last three values of the x-axis and
# y-axis in a more friendly way. It is not necessary to get the points
# themselves, it is just a presentation artifice for the article
with np.printoptions(threshold=10, precision=3, suppress=True):
display(np.array(p[0].get_points()))
array([[ 0.1 , 0.164, 0.234, ..., 3.836, 3.921, 4. ],
[244.654, 149.528, 104.515, ..., 6.378, 6.24 , 6.116]])
As we can see, there are the values of the “x axis” from 0.1 to 4 in the first list, as requested when calling the plot method, and the associated y values in the second list. Let’s see the number of elements in each list:
len(p[0].get_points()[0]), len(p[0].get_points()[1])
\(\displaystyle \left( 65, \ 65\right)\)
The number of points is estimated by SymPy itself. It is possible to customize, but honestly, I have never found it necessary, the quantities calculated by SymPy have always provided me with well-presented graphs.
Let’s understand the structure of the object generated by get_points. We have
already seen that there are two lists inside the object, but let’s see the type
of the object:
type(p[0].get_points())
tuple
Thus, we see that it is a tuple and, inside the tuple, there is a list for the x values and another for the y values. We already know everything we need. Let’s associate a variable with this object, create a Matplotlib figure and make our graph:
data = p[0].get_points()
fig, ax = plt.subplots(figsize=(8, 6), facecolor=(1, 1, 1))
ax.plot(*data) # unpacking the data
[<matplotlib.lines.Line2D at 0x7f7cf5aa4800>]
OK, this is what we alreary had so far… so what?
So now we can automate the whole creation process. Let’s generate a graph with 6
temperatures. If we were to do it via SymPy, it would be a tedious procedure,
creating each expression and passing them one by one to the plot method of
SymPy. Now we can do this more quickly, as we know how to extract data from
SymPy.
Let’s start by creating our temperature values. No one deserves to do math in Kelvin in your head, so let’s think of round values in Celsius and let Python do the conversion for us:
temperatures_celsius = (-100, 0, 25, 100, 500, 1000)
temperatures_kelvin = tuple(t + 273.15 for t in temperatures_celsius)
Now, think. We want six graphs, we need the data for each one. Let’s create an
empty list to store this data. And this data needs to be extracted from SymPy
expressions with get_points. We can create each expression and extract, but we
don’t actually need to create an expression for each temperature, we can reuse
the same variable, there is no interest in objects for each expression.
Basically, we describe a for loop:
plots = []
for temperature in temperatures_kelvin:
expr = volume.subs({R: gas_constant, n: 1, T: temperature})
p = plot(expr, (P, 0.1, 4), show=False)
data = p[0].get_points()
plots.append(data)
Note that inside the loop we did the same procedure described earlier in the article. But we were extracting the data at each step and appending it to our list. Thus, our list has 6 sets of data, a set of points (x, y), or better (P, V), for each temperature:
len(plots)
\(\displaystyle 6\)
Now it’s time to use the power of Matplotib. Following the composition we discussed earlier, let’s create a figure and associate axis systems with this figure, each system with one of the curves. And, of course, identify the figure and each axis, as well as put a legend. Full service:
fig, ax = plt.subplots(figsize=(12, 8), facecolor=(1, 1, 1))
for i, p in enumerate(plots):
ax.plot(*p, label=f'T = {temperatures_celsius[i]} °C', linewidth=3)
ax.set_ylim(0, 80)
ax.legend(fontsize=14)
ax.set_ylabel('Volume / liter', fontsize=16)
ax.set_xlabel('Pressure / atm', fontsize=16)
ax.set_title("Molar volume - ideal gas", fontsize=18)
plt.show()
That’s it! The procedure can now be extracted to a script, for example, and reproduced for any other case, not necessarily with the ideal gas law but with any other equation with only one variable.
But 2D is for the weak, we want 3D!
3D plots #
Let’s remember the objects created at the beginning:
ideal_gas_law
\(\displaystyle P V = R T n\)
volume
\(\displaystyle \frac{R T n}{P}\)
At the beginning, we fixed the temperatures at some values and generated 2D graphs, which are isotherms. But now we want the temperature to be one of the axes of the graph. Thus, the only variable we will fix is the amount of matter in moles, at the value of 1 mol:
molar_volume = volume.subs({R: gas_constant, n: 1})
molar_volume
\(\displaystyle \frac{0.08205736608096 T}{P}\)
At the beginning, we have imported the plot3d method from SymPy. The way to
use it is very similar to what we saw for plot, but now we will pass intervals
not only for pressure but also for temperature:
plot3d(molar_volume, (P, 0.1, 4), (T, 0, 1000))
<sympy.plotting.plot.Plot at 0x7f7cf595e000>
This way, we generated a PVT surface with pressure in the range between 0.1 and 4 atm and temperature in the range of 0 to 1000 K. All the two-dimensional graphs you have seen when studying gases are nothing more than cuts in planes of surfaces like this, keeping something constant: temperature (isotherm graph), pressure (isobar graph) or volume (isometric graph).
Note that SymPy already seeks to give some touches to the graph, for example, using a color gradient on the surface, indicating with warmer colors higher values of the z axis, in our case the volume axis. But we want to have even more customization power, so let’s extract the data. The procedure is analogous to what we saw earlier:
p = plot3d(molar_volume, (P, 0.1, 4), (T, 0, 1000), show=False)
p
<sympy.plotting.plot.Plot at 0x7f7cf5880fe0>
p[0]
<sympy.plotting.plot.SurfaceOver2DRangeSeries at 0x7f7cf5880b90>
Note that the object is distinct, it is SurfaceOver2DRangeSeries, literally a
surface over a range of values in two dimensions (P and T). It makes total sense
in our context. Let’s see the public attributes and methods associated with this
object:
[x for x in dir(p[0]) if not x.startswith('_')]
['end_x',
'end_y',
'expr',
'get_color_array',
'get_meshes',
'is_2Dline',
'is_3D',
'is_3Dline',
'is_3Dsurface',
'is_contour',
'is_implicit',
'is_line',
'is_parametric',
'nb_of_points_x',
'nb_of_points_y',
'start_x',
'start_y',
'surface_color',
'var_x',
'var_y']
There is no get_points and it makes sense not to have one. When dealing with
3-dimensional surface graphs, what we have are pairs of coordinates (x, y)
associated with values in z. These pairs build what is usually called a mesh
grid. Indeed, there is a method in NumPy to create
meshgrids
that you will see in several examples in the Matplotlib documentation on the
generation of 3D
surfaces.
Thus, we recognize that what we want is provided by get_meshes:
# the number of values is quite long, so I used a little bit of NumPy here
# just to show the first three and the last three values of the x-axis and
# y-axis in a more friendly way. It is not necessary to get the points
# themselves, it is just a presentation artifice for the article
with np.printoptions(threshold=10, precision=3, suppress=True, linewidth=100):
display(np.array(p[0].get_meshes()))
array([[[ 0.1 , 0.18 , 0.259, ..., 3.841, 3.92 , 4. ],
[ 0.1 , 0.18 , 0.259, ..., 3.841, 3.92 , 4. ],
[ 0.1 , 0.18 , 0.259, ..., 3.841, 3.92 , 4. ],
...,
[ 0.1 , 0.18 , 0.259, ..., 3.841, 3.92 , 4. ],
[ 0.1 , 0.18 , 0.259, ..., 3.841, 3.92 , 4. ],
[ 0.1 , 0.18 , 0.259, ..., 3.841, 3.92 , 4. ]],
[[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 20.408, 20.408, 20.408, ..., 20.408, 20.408, 20.408],
[ 40.816, 40.816, 40.816, ..., 40.816, 40.816, 40.816],
...,
[ 959.184, 959.184, 959.184, ..., 959.184, 959.184, 959.184],
[ 979.592, 979.592, 979.592, ..., 979.592, 979.592, 979.592],
[1000. , 1000. , 1000. , ..., 1000. , 1000. , 1000. ]],
[[ 0. , 0. , 0. , ..., 0. , 0. , 0. ],
[ 16.746, 9.325, 6.461, ..., 0.436, 0.427, 0.419],
[ 33.493, 18.649, 12.922, ..., 0.872, 0.854, 0.837],
...,
[ 787.081, 438.261, 303.677, ..., 20.493, 20.077, 19.677],
[ 803.827, 447.586, 310.138, ..., 20.929, 20.504, 20.096],
[ 820.574, 456.91 , 316.599, ..., 21.365, 20.931, 20.514]]])
Let’s assing to a variable:
data = p[0].get_meshes()
Now we have total control over what to do. For example, you didn’t understand the mesh part very well? Let’s make a graph with only points with the data:
fig = plt.figure(figsize=(12, 8), facecolor=(1, 1, 1))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(*data)
ax.set_xlabel('Pressure / atm')
ax.set_ylabel('Temperature / K')
ax.set_zlabel('Volume / l')
plt.show()
See how “lines” are created, sequential values in the pressure x temperature
plane, both lines parallel to the pressure axis and parallel to the temperature
axis. Note, in the representation presented above of the values provided by
get_meshes, how there are repeated values. This is the mesh of paired values
(P, T). And for each pair we have associated a value of volume V presented on
the vertical axis. I hope it has helped to understand better.
Particularly, the representation I find most useful is the “wireframe” one:
fig = plt.figure(figsize=(12, 8), facecolor=(1, 1, 1))
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(*data, rstride=5, cstride=5)
ax.set_xlabel('Pressure / atm')
ax.set_ylabel('Temperature / K')
ax.set_zlabel('Volume / l')
plt.show()
Do you recognize the isotherms and isobars in the graph?
Finally, we can also create surfaces:
fig = plt.figure(figsize=(12, 8), facecolor=(1, 1, 1))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(*data)
ax.set_xlabel('Pressure / atm')
ax.set_ylabel('Temperature / K')
ax.set_zlabel('Volume / l')
plt.show()
And color it with any colormap available in Matplotlib:
fig = plt.figure(figsize=(12, 8), facecolor=(1, 1, 1))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(*data, cmap='rainbow')
ax.set_xlabel('Pressure / atm')
ax.set_ylabel('Temperature / K')
ax.set_zlabel('Volume / l')
plt.show()
And, since we are here, let’s make this graph rotate. Why? Because yes, the graph is mine and I do wherever I want ;-)
%matplotlib notebook
fig = plt.figure(figsize=(6, 4), facecolor=(1, 1, 1))
ax = fig.add_subplot(111, projection="3d")
ax.plot_surface(*data, cmap="rainbow")
ax.set_xlabel("Pressure / atm")
ax.set_ylabel("Temperature / K")
ax.set_zlabel("Volume / l")
def animate(angle):
ax = plt.gca() # gca = get current axis
ax.view_init(30, angle)
anim = animation.FuncAnimation(
fig, animate, frames=361, interval=1, repeat=False # 360° rotation
)
plt.show()
<IPython.core.display.Javascript object>
Conclusion #
Two of the libraries I like the most in the same article. Although this article is part of the series on SymPy, much has been discussed about the way to build graphs with Matplotlib. This knowledge is essential to be able to extract the most of everything that is possible to do in terms of graphs. We also saw how SymPy behaves when interacting with other libraries, since not always internal instructions of SymPy are passed on to other libraries.
See you next time.