Basic Plotting with matplotlib

The Python ecosystem gives us the option of visualizing relationships between numbers via matplotlib, a plotting library (i.e., not part of core Python) which can produce quality figures. Inside the matplotlib package is the matplotlib.pyplot module, which is used to produce figures in a MATLAB-like environment.

Here a simple example code for the matplotlib.pyplot module.

import matplotlib.pyplot as plt

def plotex(cxs,cys,dxs,dys):
    plt.xlabel('$x$', fontsize=10)
    plt.ylabel('$f(x)$', fontsize=10)
    plt.plot(cxs, cys, 'r-', label='quadratic: $y=x^2$')
    plt.plot(dxs, dys, 'b--^', label='other function: $y=x^{1.8}-1/2$')
    plt.legend()
    plt.show()

x1 = [0.1*i for i in range(60)]
y1 = [x**2 for x in x1]

x2 = [i for i in range(7)]
y2 = [x**1.8 - 0.5 for x in x2]

plotex(x1, y1, x2, y2)

We start by importing matplotlib.pyplot in the (standard) way which allows us to use it below without repeated typing of unnecessary characters.

We then define a function, plotex(), that takes care of the plotting, whereas the main program simply introduces four list comprehensions and then calls our function. If you’re still a beginner, you may be wondering why we defined a Python function in this code. An important design principle in computer science goes by the name of separation of concerns (or sometimes information hiding or encapsulation): each aspect of the program should be handled separately. In our case, this means that each component of our task should be handled in a separate function.

Let’s discuss this function in more detail. Its parameters are (meant to be) four lists, namely two pairs of \(x_i\) and \(y_i\) values. The function body starts by using xlabel() and ylabel() to provide labels for the x and y axes. It then creates individual curves/sets of points by using matplotlib’s function plot(), passing in the x-axis values as the first argument and the y-axis values as the second argument. The third positional argument to plot() is the format string: this corresponds to the color and point/line type. In the first case, we used r for red and - for a solid line.

The fourth argument to plot() is a keyword argument containing the label corresponding to the curve. In the second call to plot() we pass in a different format string and label (and, obviously, different lists); observe that we used two style options in the format string: -- to denote a dashed line and ˆ to denote the points with a triangle marker. The function concludes by calling legend(), which is responsible for making the legend appear, and show(), which makes the plot actually appear on our screen.

We could fine-tune almost all aspects of our plots, including basic things like line width, font size, and so on. For example, we could get LaTeX-like equations by putting dollar signs inside our string, e.g., ‘$x i$’ appears as \(x_i\).

Electric Field of a Distribution of Point Charges

Very briefly, let us recall Coulomb’s law: the force on a test charge Q located at point P (at the position r), coming from a single point charge q₀ located at r₀ is given by:

\(\mathbf{F}_{0}=k \frac{q_{0} Q}{\left(\mathbf{r}-\mathbf{r}_{0}\right)^{2}} \frac{\mathbf{r}-\mathbf{r}_{0}}{\left|\mathbf{r}-\mathbf{r}_{0}\right|}\)

where Coulomb’s constant is \(k=1 /\left(4 \pi \epsilon_{0}\right)\) in SI units (and \(\epsilon_0\) is the permittivity of free space).

The force is proportional to the product of the two charges, inversely proportional to the square of the distance between the two charges, and points along the line from charge q₀ to charge Q. The electric field is then the ratio of the force F₀ with the test charge Q in the limit where the magnitude of the test charge goes to zero. In practice, this given us:

\(\mathbf{E}_{0}(\mathbf{r})=k q_{0} \frac{\mathbf{r}-\mathbf{r}_{0}}{\left|\mathbf{r}-\mathbf{r}_{0}\right|^{3}}\)

where we cancelled out the Q and also took the opportunity to combine the two denominators. This is the electric field at the location r due to the point charge q₀ at r₀.

If we were faced with more than one point charge, we could apply the principle of superposition: the total force on Q is made up of the vector sum of the individual forces acting on Q. As a result, if we were dealing with the n point charges q₀, q₁,…, q_{n −1} located at r₀, r₁,…, r_n−1 (respectively) then the electric field at the location r is:

\(\mathbf{E}(\mathbf{r})=\sum_{i=0}^{n-1} \mathbf{E}_{i}(\mathbf{r})=\sum_{i=0}^{n-1} k q_{i} \frac{\mathbf{r}-\mathbf{r}_{i}}{\left|\mathbf{r}-\mathbf{r}_{i}\right|^{3}}\)

Note that you can consider this total electric field at any point in space, r. Note, also, that the electric field is a vector quantity: at any point in space this E has a magnitude and a direction. One way of visualizing vector fields consists of drawing field lines, namely imaginary curves that help us keep track of the direction of the field. More specifically, the tangent of a field line at a given point gives us the direction of the electric field at that point. Field lines do not cross; they start at positive charges (“sources”) and end at negative charges (“sinks”).

We will plot the electric field lines in Python; while more sophisticated ways of visualizing a vector field exist (e.g., line integral convolution), what we describe below should be enough to give you a qualitative feel for things.

We are faced with two tasks: first, we need to produce the electric field (vector) at several points near the charges and, second, we need to plot the field lines in such a way that we can physically interpret what is happening.

Code below is a Python implementation, where Coulomb’s constant is divided out for simplicity.

import numpy as np
import matplotlib.pyplot as plt
from math import sqrt
from copy import deepcopy

def makefield(xs, ys):
    qtopos = {1: (-1,0), -1: (1,0)}
    n = len(xs)
    Exs = [[0. for k in range(n)] for j in range(n)]
    Eys = deepcopy(Exs)
    for j,x in enumerate(xs):
        for k,y in enumerate(ys):
            for q,pos in qtopos.items():
                posx, posy = pos
                R = sqrt((x - posx)**2 + (y - posy)**2)
                Exs[k][j] += q*(x - posx)/R**3
                Eys[k][j] += q*(y - posy)/R**3
    return Exs, Eys

def plotfield(boxl,n):
    xs = [-boxl + i*2*boxl/(n-1) for i in range(n)]
    ys = xs[:]
    Exs, Eys = makefield(xs, ys)
    xs=np.array(xs); ys=np.array(ys)
    Exs=np.array(Exs); Eys=np.array(Eys)
    plt.streamplot(xs, ys, Exs, Eys, density=1.5, color='m')
    plt.xlabel('$x$')
    plt.ylabel('$y$')
    plt.show()

plotfield(2.,20)

We start by importing numpy and matplotlib, since the heavy lifting will be done by the function streamplot(), which expects NumPy arrays as input. We also import the square root and the deepcopy() function, which can create a distinct list-of-lists.

The function makefield() takes in two lists, xs and ys, corresponding to the coordinates at which we wish to evaluate the electric field (x and y together make up r). We also need some way of storing the r_i at which the point charges are located. We have opted to store these in a dictionary, which maps from charge q_i to position r_i. For each position r we need to evaluate E(r): in two dimensions, this is made up of E_x (r) and E_y (r), namely the two Cartesian components of the total electric field. Focusing on only one of these for the moment, say E_x (r), we realize that we need to store its value for any possible r, i.e., for any possible x and y values. We decide to use a list-of-lists, produced by a nested list comprehension. We then create another list-of-lists, for Ey(r). We need to map out (i.e., store) the value of the x and y components of the total electric field, at all the desired values of the vector r, namely, on a two-dimensional grid made up of xs and ys. This entails computing the electric field (contribution from a given point charge q_i ) at all possible y’s for a given x, and then iterating over all possible x’s. We also need to iterate over our point charges q_i and their locations r_i; we do this by saying for q, pos in qtopos.items(): at which point we unpack pos into posx and posy. We thus end up with three nested loops: one over possible x values, one over possible y values, and one over i. All three of these are written idiomatically, employing items() and enumerate().

Our second function, plotfield(), is where we build our two-dimensional grid for the xs and ys. We take in as parameters the length L and the number of points n we wish to use in each dimension and create our xs using a list comprehension; all we’re doing is picking x’s from −L to L. We then create a copy of xs and name it ys. After this,
we call our very own makefield() to produce the two lists-of-lists containing E_x (r) and E_y (r) for many different choices of r.

The result of running this code is shown in the following figure.