gnuplot for dummies…

At first glance, it might seem that Gnuplot is more difficult to use for plotting than MS Excel. But it only seems that way, the entry threshold is a little higher (you can’t use a mouse to click it, you need to read the documentation here), but in practice, it comes out much easier and more convenient. I wrote a script once and have been using it all my life. It’s really much more difficult for me to plot a graph in Excel, where everything is not logical than in Gnuplot. And the main advantage of Gnuplot is that you can embed it into your programs and visualize data on the fly. Also Gnuplot without any problems builds a graph from a 30-gigabyte file of statistical data, while Exel simply crashed and could not open it.

The pluses of Gnuplot include the fact that it can be easily integrated into the code in standard programming languages. There are ready-made libraries for many languages, I personally came across Fortran and Python. Thus, you can generate graphs directly from your program.

Gnuplot — real application

There are two ways to work with Gnuplot: command mode and script execution mode. I recommend using the second mode immediately, as it is the most convenient and fastest. Moreover, the functionality is absolutely the same.

The graph is plotted by the plot command. As command parameters, you can specify a function or the name of the data file. As well as which column of data to use and how to connect the points, how to denote them, etc. Let’s illustrate.

gnuplot> plot sin(x)
Gnuplot interactive plot of sin(x) function

After executing the command, a window will open where there will be a sine graph with default signatures. Let’s improve this graph and figure out some additional options at the same time.
Let’s sign the axes.

set xlabel "X"

Specifies the label for the abscissa axis.

set ylabel "Y"

Specifies the label for the ordinate axis.

Let’s add a grid to show where the graph is built.

set grid

We will set the limits of the values to which the graph will be limited via the command

set yrange [-1.1:1.1]

Thus, the graph will be drawn from the minimum value of -1.1 to a maximum value of 1.1.

Likewise, I set the range for the abscissa so that only one period of the sinusoid is visible.

set xrange[-pi:pi]

We need to add a title to our schedule:

set title "Gnuplot's examples" font "Helvetica Bold, 20"

Note that you can set fonts and their sizes. See the Gnuplot documentation for what fonts can be used.

And finally, in addition to the sine on the graph, let’s also draw the cosine, and also set the line type and its color. And also add legends, what are we drawing.

plot sin(x) title "sin(x)" lc rgb "red", cos(x) title "cos(x)" lc rgb "green"

Here we draw two graphs on one canvas, in red and green. In general, there are a lot of options for lines (dotted line, stroke, solid), line widths, colors. Like the types of points. My goal is only to demonstrate the range of possibilities. Let’s bring all the commands into one pile and execute them sequentially.

the result:

If you honestly repeated all this after me, you might have noticed that manually entering it every time, even copying, is somehow not comme il faut. But is this a ready-made script? So let’s do it!

a gnuplot script

We make it executable and run it.

chmod +x gnuplottest.gpi

As a result, we get the same window with charts. If you don’t add -persist to the title, the window will automatically close after the script is executed.

In order to output data to a file and not to the screen, you need to reassign the terminal.

set terminal png size 800, 600
set output "result.png"
rep

As you might guess, we indicate the file type, its resolution, then we indicate the file name. Add these two lines to the beginning of our script, and we get this picture in the current folder.

Our final Gnuplot graph

In order to save to postscript, you need to use the following commands:

set terminal postscript eps enhanced color solid
set output "result.ps"
rep

In Python?

Here a simple code for plot a cosine function in Python using matplotlib:

import numpy as np
import matplotlib.pyplot as plt
t = np.arange(0.0,2.0,0.01)
s = 1+ np.cos(2*np.pi*t)
plt.plot(t,s,'*')
plt.xlabel('Time(t)')
plt.ylabel('Voltage(mV)')
plt.title('Cosine wave plot(cos(x))')
plt.show()

PyFiglet…

FIGlet – as you can read here – is a program for making large letters out of ordinary text.

A very quick note on how to add FontArt to Python’s output view.I will give a few tips on how to add FontArt to the python output display.

The method is quite easy, the first thing we have to do is install the PyFiglet library. PyFiglet converts ASCII text to ASCII art font.

pip install pyfiglet

After the installation is complete, then open the code editor and type the following code:

# You can get the type of font itself at the following link.
# http://www.figlet.org/fontdb.cgi

import pyfiglet
result = pyfiglet.figlet_format('wow', font='banner3')
print(result)

The result of the code is like this:

In addition, we can also change the displayed font, by adding a font parameter to the .figlet_format function:

# You can get the type of font itself at the following link.
# http://www.figlet.org/fontdb.cgi

import pyfiglet
result = pyfiglet.figlet_format('wow', font='colossal')
print(result)

The result of the code is like this:

Below, a complete code for the project:

import pyfiglet
import argparse

# HOW TO USE
# python asciart.py
# python asciart.py --text "Ciao"
# python asciart.py --text "Ciao" --font banner3-d

# You can see the font list here:
# http://www.figlet.org/fontdb.cgi

parser = argparse.ArgumentParser()
parser.add_argument('-t', '--text', help='text')
parser.add_argument('-f', '--font', help='font')

text = 'This is sample text'

args = parser.parse_args()
if args.text:
    text = args.text

if args.font:
    result = pyfiglet.figlet_format(text, font=args.font)
else:
    result = pyfiglet.figlet_format(text)

print(result)

The Fibonacci sequence…

Fibonacci

The Fibonacci sequence is one of the most well known mathematical sequences and is the most basic example of recurrence relations. Each number in the sequence consists of the sum of the previous two numbers and the starting two numbers are \(0\) and \(1\). It goes something like this:

\(1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144 \text{ and so on forever…}\)

We will go through a few different approaches for the optimal solution:

  1. Using simple recursion
  2. Using cache with recursion
  3. Using Binet’s formula

Using simple recursion

This is a very simple and easy way to return the nth Fibonacci number in Python:

def recursiveFib(n):
    if n == 1 or n == 2:
        return 1

    return recursiveFib(n - 1) + recursiveFib(n - 2)

It uses recursion to repeatedly call itself multiple times calculating the previous number and using it to work out the next. But that is also its downside, as the function extremely inefficient and resource intensive, this is because at every stage it calculates the previous two numbers, and the previous two numbers of those number etc. Soon you reach a point where it takes too long to calculate the next number – for example on my computer it took me 1.61 seconds to calculate the 35th number. This will obviously be extremely slow, and basically impossible, to calculate higher values.

processing time of the recoursiveFibo routine

Using cache with recursion

Since we constantly calculate the previous two numbers, we can harness the power of caching to store the number, so we do not need to calculate them anymore. The built-in functools module allows us to use a least recently used cache; a type of cache which organises the items in order of use. This can speed up the process by a huge amount.

# procedura ricorsiva con cache
@lru_cache()
def recursiveFibCached(n):
    if n == 1 or n == 2:
        return 1

    return recursiveFibCached(n - 1) + recursiveFibCached (n - 2)

Firstly, we need to import the lru_cache decorator from the functools module and place it before our function. We can supply a maxsize value to tell the cache how many items to store, but by default it is 128 and that works perfectly fine.

Using cache with recursion

Using Binet’s formula

Binet’s formula is a formula which can be used to calculate the nth term of the Fibonacci sequence, which is exactly what we want to do, and is named after it was derived by the french mathematician Jacques Philippe Marie Binet. The formula is shown below:

\(F_n = \frac{\phi^n-(-\phi)^{-n}}{\sqrt{5}}\)

where:

\(\phi = \frac{1+\sqrt{5}}{2}\)

We can convert thing formula into Python and start using it straight away:

def formulaFib(n):
    root_5 = 5 ** 0.5
    phi = ((1 + root_5) / 2)

    a = ((phi ** n) - ((-phi) ** (-n))) / root_5

    return round(a)
Using Binet’s formula

Note for the Python implementation we need to return the rounded version of the number we calculate, this is because when we calculate a large number Python will return to us a number where there can be 20+ trailing 9’s after the decimal place.

This is all well and good since now we do not have any loops and can instantly compute the answer, right? Well, there is a small catch to this method. If we try to calculate anything above the 1475th number, we will run into an error; OverflowError: (34, 'Numerical result out of range'). This is due to the way floats are implemented in Python and they can only have a certain maximum value, which we exceed using this method.

However, a fix to this is very easy to implement. We can use a built-in module called decimal to create a decimal object with a much higher precision and size to use in our equation:

import decimal

def formulaFibWithDecimal(n):
    decimal.getcontext().prec = 10000

    root_5 = decimal.Decimal(5).sqrt()
    phi = ((1 + root_5) / 2)

    a = ((phi ** n) - ((-phi) ** (-n))) / root_5

    return round(a)

In this new function we set the precision value to 10’000 digits long and we convert our root 5 value into a decimal object value and use that in our equation. This allows us to easily calculate the 10’000th number in the sequence in an astonishing 41.9 seconds, a huge improvement over all our previous methods.

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)
A matplotlib.pyplot module example

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.

Some of the most important line styles/markers in matplotlib

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 q0 located at r0 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 q0 to charge Q. The electric field is then the ratio of the force F0 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 q0 at r0.

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 q0, q1,…, qn −1 located at r0, r1,…, rn−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 ri at which the point charges are located. We have opted to store these in a dictionary, which maps from charge qi to position ri. For each position r we need to evaluate E(r): in two dimensions, this is made up of Ex (r) and Ey (r), namely the two Cartesian components of the total electric field. Focusing on only one of these for the moment, say Ex (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 qi ) 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 qi and their locations ri; 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 Ex (r) and Ey (r) for many different choices of r.

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

Visualizing the electric fields resulting from two point charges.