Plot della funzione Gamma con matplotlib

Scarica l’articolo in formato pdf

Abstract

Utilizzando alcune librerie Python di calcolo numerico utilizzeremo il pacchetto matplotlib per disegnare il grafico della funzione Gamma di Eulero. Questa infatti è la funzione scelta per introdurre questo interessante ambiente di computazione e visualizzazione, integrato e circoscritto all’interno del linguaggio Python.

Lo scopo è quello di mostrare le capacità di questo software libero e dare un idea delle caratteristiche di un linguaggio testuale basato su Python per il calcolo e la visualizzazione tecnico-scientifica.

Numeri, numeri ed ancora numeri per Python

Il linguaggio Python è particolarmente portato per i numeri. Lo dimostra questo piccolo e classico esempio in cui si calcola il fattoriale di un numero intero elevato.
Il fattoriale non è altro che il prodotto di tutti i numeri da 1 ad n e può essere calcolato da una funzione ricorsiva. Ecco l’esempio di codice:

def fact(n):
    if n==1:
        return 1
    return n * fact(n-1)

print(fact(100))

Salvato il file contenente il codice con il nome fact.py , per eseguirlo aprire una finestra di terminale e digitare $ python fact.py. Sarà immediatamente mostrato a video un numero di 158 cifre. Nello screenshot seguente vi ho riportato l’esecuzione di fact.py posizionato nel mio caso nella cartella Scrivania:

Esecuzione dello script Python per il fattoriale di 100

È preferibile riscrivere il codice senza copiarlo per rendersene meglio conto, inoltre ricordatevi che se non vogliamo digitare il percorso completo del nome del file da eseguire, la directory di lavoro del terminale o della console di Windows, deve essere la stessa di quella in cui si trova il file. Per settare la directory di lavoro, e quindi spostarsi nell’albero del file system, basta usare il comando cd (che sta per change directory) seguito dal nome della directory.

Ultima cosa da ricordare, questa volta agli utenti Linux, è l’assegnazione agli script Python dei permessi di esecuzione, per esempio semplicemente spuntando il flag Esecuzione nella scheda Permessi del dialogo Proprietà del file stesso.

Preparazione del sistema

In Ubuntu installare i file necessari per svolgere l’esercizio proposto è semplice perché l’interprete Python è pre-installato di default mentre le librerie matplotlib e scipy si installano sul sistema con i due comandi (la libreria numerica numpy è una dipendenza obbligatoria del pacchetto python-matplotlib e non è necessario richiederne l’installazione esplicitamente):

$ sudo apt-get install python-matplotlib
$ sudo apt-get install python-scipy

Per gli utenti Windows serve un po’ più di tempo ma mutatis mutandis, non è difficile seguire le procedure d’installazione di Python e delle librerie numeriche indicate. La documentazione è ottima.

matplotlib plotting example

La libreria matplotlib funziona allo stesso modo di Matlab: la curva di una funzione è disegnata unendo un insieme di punti del piano le cui coordinate sono definite da un set di valori del dominio della funzione e dai corrispondenti valori calcolati nel codominio.

Quello che dobbiamo fare è scrivere un programma in Python che utilizzi motori numerici e di visualizzazione sotto forma di librerie

Consideriamo un esempio di plot di funzione per . Il set di valori sull’asse x saranno i 100 punti egualmente intervallati da 0 a : setx = np.linspace(0, 2*np.pi, 100), dove np è l’oggetto che fa riferimento alla libreria numpy definito con un precedente import numpy as np (il nome assegnato è una sorta di scorciatoia).
Il plottaggio avviene utilizzando il metodo plot dell’oggetto matplotlib.pyplot, a cui passiamo il set di valori sull’asse x ed i corrispondenti ordinate sull’asse y calcolate con la funzione seno al quadrato: plt.plot(setx, np.sin(setx)**2) dove plt è l’oggetto che fa riferimento alla libreria matplotlib.pyplot definito con import matplotlib.pyplot as plt.

Non rimane che esercitarci sul codice completo del programma che genererà un file pdf, quindi un file di qualità vettoriale che scalato non perderà risoluzione:

#!/usr/bin/python

import numpy as np
import matplotlib.pyplot as plt

setx = np.linspace(0.0, 2*np.pi, 100)
plt.plot(setx, np.sin(setx)**2)

plt.savefig('sen2plt.pdf')

Il grafico della funzione seno al quadrato (click to download the pdf)

Sul manuale troverete le usuali proprietà per personalizzare il grafico come etichette compreso testo in codice LaTeX per il typesetting matematico, spessori e colori, griglie, ecc. Le lavorazioni sui sorgenti assomigliano operativamente a quelle degli altri tool di visualizzazione come Matlab stesso o Gnuplot, ma con i vantaggi di un linguaggio ad oggetti intuitivo. Occorre solo l’abitudine ad orientarsi tra i moduli delle librerie.

Esiste anche la possibilità di visualizzare il grafico in una finestra grafica. Non seguiremo questa strada perché interattiva per l’utente e quindi poco adatta ad essere descritta in una guida. Dico solo che per usufruire di questa modalità basta sostituire nell’ultima istruzione il metodo savefig(‘nomefile’) con show().

The Gamma function

La funzione fu trovata da Eulero su sollecitazione di Christian Goldbach e rappresenta il fattoriale di un numero reale (od anche complesso). Si tratta quindi di una sorta di estensione del calcolo del fattoriale di un numero intero infatti per essi sussiste la seguente relazione:

La definizione standard della funzione Gamma di Eulero è la seguente:

Per conoscere ulteriori dettagli matematici visitate la classica pagina di Wikipedia o quella di MathWorld.

matplotlib plots Gamma

La funzione è stata implementata nel modulo speciale della libreria scipy. Il nostro obbiettivo meno ambizioso quanto a visualizzazioni 3D nel piano complesso, è tracciare il grafico per un dominio reale ed evidenziare i punti sulla curva corrispondenti ai fattoriali dei numeri interi.
Ma la funzione non si domina facilmente. Giustamente matplotlib non disegna punti a distanza infinita, quindi occorre gestire le discontinuità ovvero i valori infiniti che essa assume per 0 e per i numeri interi negativi, appiattendo i valori della funzione quando superano un range scelto opportunamente e molto al di fuori dei limiti del grafico (questo per la verità va fatto per qualsiasi altra funzione con asintoti verticali nell’intorno dei quali i valori diventano ovviamente elevatissimi).
Poi occorre infittire molto i punti in cui calcolare altrimenti nel grafico compariranno dei rami parassiti nei pressi degli asintoti verticali.

Passo infine la parola al codice:

#!/usr/bin/python

# import delle librerie necessarie
import numpy as np
from scipy.special import gamma as Gamma
import matplotlib.pyplot as plt

# definizione dei limiti del grafico
# per il facile 'aggiustamento'
x0, x1 = -5.0, 5.2
y0, y1 = -6.5, 7.0

# dominio lineare sull'asse x
x = np.linspace(x0, x1, 50000)

# dominio numeri interi
xx = [2,3,4]

# creazione dei vettori di ordinata
y = Gamma(x)
yy = Gamma(xx)

# 'gestione' delle discontinuita
ulim = 100*y1
llim = 100*y0
y[y>ulim] =  np.inf
y[y<llim] = -np.inf
# plottaggio della Gamma
plt.plot(x, y)

# aggiunta dei punti notevoli del fattoriali
# l'opzione 'ro' si riferisce a singoli caratteri:
# r = red , o = pallino
plt.plot(xx, yy , 'ro')

# limiti sugli assi
plt.axis([x0, x1, y0, y1])

# tracciamento di una griglia base
plt.grid()

# creazione del file pdf con il grafico
plt.savefig('gammaplt.pdf')

Grafico della Gamma di Eulero (click to download the pdf file)

Conclusioni

Secondo me bene si fa a tener conto delle librerie Python introdotte in questo post. Basta dare un occhiata ai numerosi talk ai meeting europei di Scipy.
Impara l’arte e mettila da parte…