Calcolo Numerico in Python

Calcolo Numerico in Python

Fabio Marcuzzi, Maria Rosaria Russo

Dipartimento di Matematica Pura ed ApplicataUniversita‘ di Padova

Italye–mail: [email protected], [email protected]

Versione 3.0

6 Giugno 2005

F.Marcuzzi - M.R.Russo

Contents

1 Introduzione 4

1.1 Sviluppo di programmi in Python . . . . . . . . . . . . . . . . . . . . . . . 4

1.2 Documentazione su Python . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 Moduli aggiuntivi per il calcolo scientifico . . . . . . . . . . . . . . . . . . 5

2 Modalita di esecuzione dei programmi scritti in Python 6

2.1 Caricamento di Moduli aggiuntivi (import e reload) . . . . . . . . . . . . 6

2.2 Lavorare con la shell interattiva di Python . . . . . . . . . . . . . . . . . . 7

2.2.1 Ripetizione dei comandi . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2.2 Interazione con le directories . . . . . . . . . . . . . . . . . . . . . . 7

2.2.3 Help in linea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3 Presentazione sintetica del Linguaggio 8

3.1 Tipi di dato . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.2 Istruzioni composte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3.2.1 La decisione if-then-else . . . . . . . . . . . . . . . . . . . . . . . . 10

3.2.2 Il ciclo for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3.2.3 Il ciclo while . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.3 Funzioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

3.4 I/O da file e da schermo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3.4.1 I/O da file . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3.4.2 I/O da schermo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.5 misura dei tempi di esecuzione . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Operazioni con Vettori e Matrici 14

4.1 Basic Linear Algebra Subroutines (BLAS) e progetto ATLAS . . . . . . . . 16

4.1.1 Operazioni di tipo BLAS 1 . . . . . . . . . . . . . . . . . . . . . . . 16



4.2 Matrici Sparse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4.2.1 Operazioni implementate in modo dedicato per il formato LL . . . . 18

4.2.2 Operazioni implementate in modo ottimizzato per il formato CSR

ed SSS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2


5 Creazione di Grafici 19

5.1 Grafici bi-dimensionali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

5.2 Grafici tri-dimensionali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5.3 Immagini . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

6 Libreria di Algebra Lineare Numerica 22

6.1 Interfaccia di Numeric alla LAPACK . . . . . . . . . . . . . . . . . . . . . 22

6.2 Soluzione di sistemi lineari con matrice sparsa . . . . . . . . . . . . . . . . 22

7 Libreria per la soluzione di Equazioni Nonlineari 22

8 Libreria per l’Interpolazione ed Approssimazione di funzioni e di dati 22

8.1 Interpolazione polinomiale . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

8.2 Approssimazione ai Minimi Quadrati Discreti . . . . . . . . . . . . . . . . 23

9 Libreria per l’Integrazione Numerica 24

10 Libreria per la soluzione di Problemi ai Valori Iniziali 24

11 Libreria per la soluzione di Problemi ai Limiti 24

12 Libreria per la Geometria Computazionale 24

13 Programmazione ad Oggetti in linguaggio Python 24

13.1 Classi e Oggetti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

13.1.1 Esempio: la classe Matrice . . . . . . . . . . . . . . . . . . . . . . . 25

3


1 Introduzione

Python e un linguaggio di scripting ad oggetti. Come proprieta fondamentale, consente

la ”tipizzazione dinamica” (dynamic typing), cioe:

1. riconosce automaticamente oggetti quali numeri, stringhe, liste, ..., e quindi non

richiede di dichiararne il tipo e la dimensione prima dell’utilizzo;

2. effettua in modo automatico l’allocazione e la gestione della memoria.

Questa caratteristica, presente in ambienti come il Matlab, contribuisce in modo

sostanziale a velocizzare la prototipazione di algoritmi di calcolo numerico.

Python e disponibile al sito www.python.it (www.python.org e il sito internazionale).

La versione a cui facciamo qui riferimento e dalla 2.3.4 in poi (in questo momento la piu

recente e la 2.4.1 ).

1.1 Sviluppo di programmi in Python

L’esecuzione di programmi scritti in Python avviene tramite la sua shell interattiva, che

viene semplicemente invocata con il comando python, sia in Linux che in DOS/Windows.

In generale, e utile dotarsi di un’IDE (Integrated Development Environment, ovvero

un ambiente di sviluppo integrato). A questo proposito ci sono varie soluzioni. Un’IDE

disponibile per Linux e Windows e Eclipse, per la quale e disponibile un plugin per

Python. Piu semplicemente:

• per Linux e conveniente utilizzare IDLE, disponibile al sito

http://www.python.org/idle/ .

• per Windows e conveniente utilizzare l’IDE ”Pythonwin”, disponibile al sito

http://starship.python.net/crew/mhammond/win32/ dentro alle Python for Win-

dows Extensions.

1.2 Documentazione su Python

Documentazione relativa a Python:

• un tutorial per iniziare e disponibile al sito

http://www.python.it/doc/Easytut/easytut-it/index.html .

• al sito http://www.python.it/ si trova una notevole quantita di informazioni in

generale sul linguaggio Python.

4


1.3 Moduli aggiuntivi per il calcolo scientifico

L’utilizzo di Python per l’implementazione di algoritmi numerici e possibile principal-

mente grazie alle sue estensioni, tra cui segnaliamo, in particolare:

• Numeric (ci riferiamo alla versione 23.8 disponibile al sito

http://sourceforge.net/project/showfiles.php?group id=1369&package id=1351

): e l’estensione che offre gli strumenti per lavorare con vettori e ma-

trici, di importanza fondamentale; documentazione ed informazioni

sono disponibili al sito http://numeric.scipy.org/ , in particolare

http://numeric.scipy.org/numpydoc/numdoc.htm .

• SciPy ( ci riferiamo alla versione disponibile al sito http://www.scipy.org/ ): e

un’estensione che offre vari strumenti di calcolo scientifico; si appoggia a Numeric.

• PySparse ( ci riferiamo alla versione disponibile al sito

http://pysparse.sourceforge.net/ ): e l’estensione che offre l’implementazione

per Python delle matrici sparse, e di alcuni algoritmi per la risoluzione di sistemi

lineari sparsi; si appoggia a Numeric.

• Pygist ( ci riferiamo alla versione disponibile al sito http://bonsai.ims.u-

tokyo.ac.jp/%7Emdehoon/software/python/pygist.html ): e l’estensione

che offre gli strumenti per creare grafici 2D (ed anche 3D);

il manuale online e disponibile al sito http://bonsai.ims.u-

tokyo.ac.jp/%7Emdehoon/software/python/pygist html/pygist.html .

• Gnuplot-py ( ci riferiamo alla versione disponibile al sito http://gnuplot-

py.sourceforge.net/ ): e un’estensione alternativa a Pygist, che interfaccia Python

al popolare Gnuplot.

• Gmsh ( ci riferiamo alla versione disponibile al sito http://www.geuz.org/gmsh/ ):

e un’estensione per la creazione di griglie 1,2,3-dimensionali irregolari (adatte ad es.

alla simulazione mediante il Metodo degli Elementi Finiti).

5


Parte I Programmazione

2 Modalita di esecuzione dei programmi scritti in Python

Ci sono due modalita principali per avviare un programma, scritto in linguaggio Python,

dal prompt della shell di Linux:

1. impartire il seguente comando (esegue il programma e ritorna alla shell di Linux):

python <nome file programma>.py

2. invocare l’ambiente basilare di sviluppo IDLE e poi lanciare il programma:

• idle

(viene aperta la finestra di IDLE)

• selezionare File–>Open e poi selezionare il file contenente il programma;

(viene aperta la finestra dell’editor con il file selezionato)

• sulla finestra dell’editor con il programma da avviare, selezionare Run–

>RunModule

(sulla shell di IDLE appaiono i risultati dell’esecuzione e rimangono disponi-

bili).

2.1 Caricamento di Moduli aggiuntivi (import e reload)

In entrambe le modalita di avvio di un programma, precedentemente citate, l’ambiente di

esecuzione viene inizializzato in modo minimale: tutti i moduli built-in e standard sono

disponibili, ma non sono stati inizializzati, con l’eccezione di sys (vari servizi di sistema),

builtin (funzioni built-in, eccezioni e None) e main .

Di conseguenza, prima di eseguire una funzione e necessario verificare di aver importato

(mediante l’istruzione import) il modulo che la contiene.

Nota bene: la chiamata della funzione varia a seconda di come si e importato il modulo.

Ad esempio, per chiamare la generica funzione <nome funzione>() del generico modulo

<nome modulo>, e possibile agire in uno dei seguenti modi:

• caricamento di un modulo mantenendone il nome:

import <nome modulo>

chiamata: <nome modulo>.<nome funzione>()

6


• caricamento di un modulo nascondendone il nome:

from <nome modulo> import *

chiamata: <nome funzione>()

• caricamento di alcune funzioni di un modulo, nascondendone il nome:

from <nome modulo> import <nome funzione> [, <altro nome funzione>, ...]

chiamata: <nome funzione>()

Per vedere i comandi disponibili nei moduli standard e possibile accedere al sito ufficiale

http://www.python.org/doc/current/lib/lib.html .

Se un modulo gia importato e stato successivamente modificato, e necessario caricarlo

nuovamente nell’ambiente. In questo caso, pero, l’istruzione import non farebbe nulla (in

quanto il modulo e gia stato importato), mentre e necessario utilizzare l’istruzione reload:

reload(<nome modulo>) .

2.2 Lavorare con la shell interattiva di Python

Dalla shell interattiva di Python e possibile impartire qualunque istruzione eseguibile

dall’interprete Python. Quindi: le istruzioni dei programmi possono essere singolarmente

provate sulla shell (ad es. per verificarne l’esecuzione). Un intero programma pero non va

scaricato direttamente sulla shell (es. mediante cut-and-paste), ma salvato su un file di

testo ed eseguito come script.

2.2.1 Ripetizione dei comandi

E’ utile ricordare che un’istruzione precedentemente impartita alla shell di Python, nella

sessione corrente, puo essere facilmente richiamata:

• con IDLE (Linux): muovendosi con le combinazioni di tasti <Alt-p> e <Alt-n>

nell’elenco storico progressivo memorizzato dall’ambiente;

• con la shell interattiva (Linux) o con pythonWin (Windows): muovendosi con le

combinazioni di tasti <ctrl-freccia su> e <ctrl-freccia giu> nell’elenco storico

progressivo memorizzato dall’ambiente .

Questa funzione e analoga a quella ottenuta in Matlab con <freccia-su> e <freccia-giu>)

.

2.2.2 Interazione con le directories

E’ necessario che il modulo os (funzioni del sistema operativo) venga prima importato:

import os

7


• Nome della directory corrente:

os.getcwd()

• Cambio della directory corrente:

os.chdir(’<nome dir figlia>’) # passa alla directory figlia specificata

os.chdir(’..’) # passa alla directory padre

os.chdir(’<full path>’) # passa alla directory specificata

• Elenco del contenuto di una directory:

os.listdir(’.’) # mostra la directory corrente

os.listdir(’..’) # mostra la directory padre

os.listdir(’<nome dir figlia>’) # mostra la directory figlia specificata

os.listdir(’<full path>’) # mostra la directory specificata

• Creazione di una directory:

os.mkdir(’<nome dir figlia>’) # crea la directory figlia specificata

os.mkdir(’<full path>’) # crea la directory specificata

2.2.3 Help in linea

Si utilizza il comando help, che fornisce informazioni testuali su:

• un modulo intero:

help(<nome modulo>)

• una funzione di un modulo:

help(<nome modulo>.<nome funzione>)

Nota: la funzione help() richiede che il modulo <nome modulo> sia gia stato importato.

3 Presentazione sintetica del Linguaggio

3.1 Tipi di dato

• string:

s = "Ciao" # creazione di una stringa

print len(s) # lunghezza in caratteri di una stringa

>>> 4

print s[0:2] # estrazione di una sotto-stringa

8


>>> ’Ci’ # NB: s[da,a] considera gli indici da "da" fino ad "a-1"

# ripetizione di una stringa:

s = "Ciao"

print s * 3

>>> ’CiaoCiaoCiao’

# concatenazione di due o piu’ stringhe:

s1 = "Ciao"

s2 = "Mondo"

print s1 + " " + s2

>>> ’Ciao Mondo’

• int: Python supporta i numeri interi a 32 bit, rappresentati mediante la notazione

complemento a 2.

• double: Python supporta i numeri in virgola mobile in doppia precisione (la singola

precisione non e prevista).

ATTENZIONE: se si scrive una costante senza il punto (”.”) viene rappresentata

come numero intero. Ad esempio, ”1” e un int, mentre ”1.” e un double, per cui si

ha:

>>> 1/2

0

>>> 1./2

0.5

>>> 1/2.

0.5

# arrotondamenti:

y = ceil(x) # "y" e’ il ’double’ pari al piu’ piccolo intero >= x

y = floor(x) # "y" e’ il ’double’ pari al piu’ piccolo intero <= x

y = round(x) # "y" e’ il ’double’ pari all’intero piu’ vicino a "x"

• complex: Python supporta i numeri complessi in virgola mobile. Esempi:

z = complex(3.1,2.2) # crea il numero complesso "3.1 + j2.2"

z.real # restituisce la parte reale di "z"

z.imag # restituisce la parte immaginaria di "z"

• lista:

e un tipo di dato utile per immagazzinare piu di un valore in un’unica variabile.

Vediamo le principali operazioni possibili:

# creazione di una lista (in generale, una lista e’ eterogenea):

listino = [’panino’,3.70,’acqua’,0.50,’gelato’,0.99]

listino[2] # accesso all’elemento di indice 2 (NB: gli indici partono da 0)

9


listino[2] = 3 # imposta a 3 l’elemento di indice 2

del listino[2] # cancella l’elemento di indice 2

len(listino) # restituisce la lunghezza della lista

"vino" in listino # e’ true se "vino" e’ un elemento della lista

"vino" not in listino # e’ vero se "vino" non e’ un elemento della lista

list.sort() # ordina la lista

list.index("acqua") # restituisce l’indice dove viene trovata la

# prima occorrenza di "acqua"

list.append("vino") # aggiunge l’elemento "vino" alla fine della lista

3.2 Istruzioni composte

3.2.1 La decisione if-then-else

Permette di esprimere un’alternativa in base ad una condizione da testare durante

l’esecuzione del programma :

if <condizione>:

<istruzioni ....>

else:

<istruzioni ....>

#end

Dal punto di vista del calcolo numerico, questo e un costrutto poco favorevole poiche, in

generale:

• nella valutazione della condizione non vengono eseguiti propriamente dei conti;

• non favorisce la creazione di un flusso omogeneo di calcoli (condizione ideale per

ottenere una buona velocita di esecuzione nelle macchine di calcolo moderne).

Dunque, esso va messo in un algoritmo solo se strettamente necessario e preferibilmente

non all’interno di un ciclo.

3.2.2 Il ciclo for

Permette di ripetere un certo insieme di istruzioni per un numero prefissato di iterazioni

:

for <variabile> in <vettore di interi>:

<istruzioni ...>

#end

10


Il <vettore di interi> puo essere espresso mediante una lista oppure un array (Nu-

meric.array), e <variabile> assume al passo i-esimo il valore della sua i-esima componente

. Esempi:

1. for i in range(1,2000,10), che significa ” i parte da 1, si incrementa con

passo 10 ed il ciclo si arresta quando i >= 2000” (in linguaggio C, sarebbe

for (i=1; i<2000; i+=10) {...});

2. for i in [1 2 4 6 25], che significa ”il ciclo esegue 5 iterazioni; alla prima i vale

1, alla seconda 2, alla terza 4, alla quarta 6, alla quinta 25”.

3. for i in stringa, che significa ” i assume valore pari ai singoli caratteri della

”stringa”, partendo dal primo e fino all’ultimo”.

Questo costrutto e largamente utilizzato nel calcolo numerico e si presta a vari tipi di

ottimizzazioni, sia su calcolatori sequenziali che paralleli.

Tipicamente la <variabile> assume valori da 1 a N , dove N e una dimensione fonda-

mentale del problema, ad es. il numero di componenti di un vettore o di righe (o colonne)

di una matrice.

3.2.3 Il ciclo while

Permette di ripetere un insieme di istruzioni per un numero di iterazioni condizionato dal

soddisfacimento di una condizione di arresto :

while <variabile> == <espressione>:

<istruzioni ...>

#end

Questo costrutto e frequentemente utilizzato nel calcolo numerico, ad es. per

l’implementazione di metodi iterativi (largamente utilizzati per la risoluzione di varie

classi di problemi numerici), e cioe di metodi che, partendo da un’approssimazione iniziale

della soluzione esatta, la migliorano per incrementi successivi utilizzando un metodo di

dimostrata convergenza.

3.3 Funzioni

Un modulo Python (lo identifichiamo qui con un file, <nome file>.py) puo contenere una

o piu funzioni al suo interno. Importando il modulo, si rendono disponibili le funzioni in

esso contenute.

In generale, una funzione la si definisce nel seguente modo:

def <nome funzione>([<parametro 1>, <parametro 2>, ...]):

<istruzioni...>

<...con eventuale "return">

11


L’identificativo def (abbreviazione di ”define”, definisci) inizializza la definizione di una

funzione. La presenza di parametri nella definizione di una funzione e opzionale, cioe si

puo definire una funzione che non riceve alcun parametro in ingresso.

Le istruzioni dopo il ”:” vengono eseguite ogni volta che la funzione viene chiamata dal

programma e sono indentate rispetto a ”def”. Quando l’indentazione torna al livello della

funzione o quando incontra l’istruzione return (che e quindi opzionale), la funzione ter-

mina.

Le variabili passate come argomenti dal programma chiamante non vengono modificate

dalla funzione (chiamata per valore).

Le funzioni permettono alle variabili locali di esistere solo ed unicamente all’interno della

funzione stessa, e di nascondere alle istruzioni di questa le altre variabili definite all’esterno

della funzione.

3.4 I/O da file e da schermo

3.4.1 I/O da file

• Scrittura di un file di testo:

out_file = open("prova.txt","w") # apertura del file in scrittura ("w")

out_file.write("Prima riga di testo\nSeconda riga di testo\n...\n")

# scrittura del file: "\n" manda a capo sulla riga successiva

out_file.close() # chiusura del file

• Lettura di un file di testo in una sola operazione:

in_file = open("prova.txt","r") # apertura del file in lettura ("r")

testo = in_file.read() # lettura del file e memorizzazione

# nella variabile string ’testo’

in_file.close() # chiusura del file

• Lettura di un file di testo riga-per-riga:

in_file = open("prova.txt","r") # apertura del file in lettura ("r")

riga = in_file.read() # lettura del file

while len(riga) > 0: # la funzione readline() ritorna una stringa vuota

# quando viene raggiunta la fine del file.

riga = in_file.read() # leggo la prossima riga del file di testo

<altre istruzioni ...>

#end

in_file.close() # chiusura del file

Nota: l’estensione del file puo essere qualunque, e non necessariamente ”txt”.

12


3.4.2 I/O da schermo

• stampa su schermo:

print "2 + 2 = ", 2+2

>>> 2 + 2 = 4

• richiesta di un numero da tastiera:

num = input("scrivi un numero:")

• richiesta di una stringa da tastiera:

str = raw_input("scrivi una stringa:")

3.5 misura dei tempi di esecuzione

E’ necessario che il modulo time venga prima importato:

import time

dopodiche l’istante attuale, misurato in secondi, viene preferibilmente fornito dalla

funzione:

• in Linux:

istante_1 = time.time() # l’istante attuale viene memorizzato

# nella variabile "istante_1"

• in Windows:

istante_1 = time.clock() # l’istante attuale viene memorizzato

# nella variabile "istante_1"

Quindi, e molto semplice misurare l’intervallo di tempo necessario all’esecuzione di un

insieme di istruzioni (al limite un intero programma): basta calcolare la differenza tra

finale ed istante iniziale.

Attenzione, pero: il programma scritto in Python non viene compilato in linguaggio

macchina (come avviene per i programmi scritti in linguaggio C/C++ o in Fortran),

bensi’ compilato in bytecode (analogamente a quanto avviene in Java) e questo poi viene

interpretato ed eseguito. Questo significa che:

13


1. il tempo di esecuzione, misurato da time.clock(), puo essere in gran parte

dovuto ad operazioni esterne all’algoritmo di calcolo numerico in esecuzione, e

che quindi prendere il tempo di esecuzione come tempo di calcolo puo risultare

un’approssimazion alquanto grossolana;

2. lo stile di programmazione incide sul tempo di esecuzione molto di piu di quanto

avvenga con i linguaggi compilati (es. C/C++ e Fortran); e necessario dunque fare

la dovuta attenzione. Ad esempio, evitare il piu possibile di fare cicli (for oppure

while) innestati.

4 Operazioni con Vettori e Matrici

L’elaborazione di vettori e matrici non puo essere fatta con gli strumenti base di Python,

ad es. liste e cicli for, in quanto risulterebbe assai lenta se confrontata con implementazioni

in linguaggio C o Fortran. Per questo motivo e stata sviluppata l’estensione Numeric, detta

anche NumPy, che consente di eseguire operazioni su vettori e matrici con un’efficienza

paragonabile a quella che si ha con i linguaggi compilati (ad es. C e Fortran, appunto).

La documentazione HTML relativa all’estensione Numeric puo essere trovata al sito:

Vediamo una sintesi delle operazioni principali:

from Numeric import * # caricamento del modulo Numeric

• Creazione di Vettori e Matrici:

v = zeros(n) # crea un vettore "v" di "n" componenti

# di tipo ’int’ e di valore zero

v = zeros(n, ’Double’) # crea un vettore "v" di "n" componenti

# di tipo ’complex’ e di valore ’0+j0’

v = zeros(n, ’double’) # crea un vettore "v" di "n" componenti

# di tipo ’double’ e di valore zero

v = ones(n, ’double’) # crea un vettore "v" di "n" componenti

# di tipo ’double’ e di valore zero

M = zeros((m,n), ’double’) # crea una matrice "M" di "m" righe ed "n" colonne,

# di tipo ’double’ e di valore uno

M = ones((m,n), ’double’) # crea una matrice "M" di "m" righe ed "n" colonne,

# di tipo ’double’ e di valore uno

v = arrayrange(a, b, delta, ’double’) # crea il vettore di double:

# (a, a+delta, ..., max(a+x*delta < b))

v = arange(a, b, delta, ’double’) # ’arange’ e’ sinonimo di ’arrayrange’

Nota: il tipo double corrisponde al numero in virgola mobile in doppia precisione

(il double del linguaggio C).

• Conversione da/in liste:

pl = [1.1, -0.9, 5, 23]; v = array(pl) # crea il vettore "v" dalla lista "pl"

r1 = [1.1, -0.9, 5]; r2 = [-3.1, 5, 2.1]; M = array([r1, r2]) # crea "M"

14


• Interazione con la dimensionalita di Vettori e Matrici:

v = array([2, 0, 3.1, 0.1, 6, 13])

v.shape # restituisce una tupla con la lunghezza di ogni dimensione di "v"

v.shape = (2,3) # trasforma "v" in una matrice 2x3

v = transpose(v) # assegna a ’v’ (che adesso e’ una matrice) la sua trasposta

size(v) # restituisce il numero totale di elementi di "v"

v.shape = (size(v),) # trasforma "v" nuovamente in un vettore di 6 elementi

len(v) # restituisce la lunghezza della prima dimensione di "v"

• Indicizzazione di Vettori e Matrici:

v = arange(-1, 0.25, 3.1) #

v[-1] = v[0] # setta l’ultimo elemento al valore del primo

v[2:5] = 0.28 # setta v[2], v[3] e v[4] al valore 0.28

M = array([[-1, 0.25, 3.1],[0.1, 0.3, -0.55]]) #

print M[0,1] # stampa il valore dell’elemento di indici 0,1

M[r,c] = 9.77 # assegna l’elemento di indici r,c al valore 9.77

print M[:,1] # stampa la seconda colonna di ’M’

M[:,:] = 1.1 # assegna tutti gli elementi di ’M’ al valore 1.1

b = M[:,2] # NB:crea un riferimento alla terza colonna di ’M’,

# non una copia !

# Quindi, modificando ’b’ si modifica anche ’M’ e viceversa

b = M[:,2].copy() # NB: crea una copia della terza colonna di ’M’

• Calcoli con Vettori e Matrici:

NB: fare cicli (for o while) sugli elementi di un vettore o di una matrice dovrebbe

essere evitato, per non rallentare troppo l’esecuzione. E’ consigliabile utilizzare il

piu possibile le molte operazioni vettoriali e matriciali implementate in Numeric,

presentate qui nel seguito.

b = v*5.5 - 1.0 # crea il vettore ’b’ in cui ogni componente e’ pari alla

# corrispondente nel vettore ’v’, moltiplicata per 5.5

# e a cui viene sottratto il valore 1.0

# operatori aritmetici vettoriali "sul posto" (analogo per le matrici):

v *= 1.155 # moltiplica per 1.155 ogni componente di ’v’

v -= 0.1 # sottrae 0.1 ad ogni componente di ’v’

v /= 1.44 # divide per 1.44 ogni componente di ’v’

v += 3.1 # somma 3.1 ad ogni componente di ’v’

v **= 2 # eleva alla potenza 2 ogni componente di ’v’

# le funzioni trigonometriche,... sono applicabili a vettori e matrici:

v = cos(v) # ogni componente di ’v’ viene sostituita dal suo coseno

b = sin(v) # crea il vettore ’b’ in cui ogni componente e’ pari al

# seno della componente corrispondente del vettore ’v’

b = log(v) # (idem)... logaritmo naturale

b = exp(v) # (idem)... esponenziale

b = sqrt(v) # (idem)... radice quadrata

15


Nota: l’estensione Numeric dispone anche di un modulo, ’MLab.py’, che implementa

una serie di funzioni con una sintassi compatibile con Matlab (cfr. la documen-

tazione di Numeric) :

from MLab import *

v = rand(n) # creazione di un vettore random di lunghezza "n"

M = rand(m,n) # creazione di una matrice random di lunghezza "mxn"

L = tril(M) # crea la matrice "L" con la parte triangolare inferiore di "M"

U = triu(M) # crea la matrice "U" con la parte triangolare superiore di "M"

• Verifica e modifica del tipo di dato:

v = array([1.2, 3.3, 2])

isinstance(v, ArrayType) # dice se ’v’ e’ un oggetto di tipo ’array’

print v.typecode() # restituisce:

# ’d’ per ’double’

# ’D’ per ’complex’

# ’i’ per ’int’

# conversione di formato:

v = v.astype(int) # converte ’v’ in formato ’int’ per troncamento

v.astype(int) # restituisce una copia di ’v’ convertita in ’int’

4.1 Basic Linear Algebra Subroutines (BLAS) e progetto ATLAS

4.1.1 Operazioni di tipo BLAS 1

• prodotto interno (o prodotto scalare) tra due vettori v (riga) e w (colonna) di

lunghezza n:

c = innerproduct(v, w)

dove c risulta essere uno scalare, pari a:

c = v[0]*w[0] + ... + v[n-1]*w[n-1]

• somma dei vettori v e w, di lunghezza n, scalati di ’alpha’ e ’beta’:

b = alpha*v + beta*w

dove b risulta essere un vettore, di lunghezza n, il cui i-esimo componente e pari a:

b[i] = alpha*v[i] + beta*w[i]

• prodotto per componenti dei vettori v e w, di lunghezza n:

b = v * w


b[i] = v[i] * w[i]

16



• prodotto esterno tra due vettori v (colonna) e w (riga), di lunghezza n:

C = outerproduct(v, w)

dove C risulta essere una matrice, di dimensione n × n, la cui componente di

posizione (i, j) e pari a:

C[i,j] = v[i]*w[j]

• prodotto matrice-vettore (matrice A di dimensione n × n e vettore colonna v di

lunghezza n o matrice di dimensione n × 1 ):

b = matrixmultiply(A, v)


b[i] = A[i,0]*v[0] + A[i,1]*v[1] + ... + A[i,n-1]*v[n-1]

• prodotto di matrici (di dimensione n × n) per componenti :

C = A * B


posizione (i, j) e pari a:

C[i,j] = A[i,j]*B[i,j]


• prodotto di matrici (di dimensione n × n):

C = matrixmultiply(A, B)


posizione (i, j) e pari al prodotto interno della riga i-esima di ’A’ per la colonna

j-esima di ’B’ :

C[i,j] = A[i,0]*B[0,j] + A[i,1]*B[1,j] + ... + A[i,n-1]*B[n-1,j]

4.2 Matrici Sparse

La matrice sparsa non e ovviamente un tipo di dato standard del linguaggio Python. La

troviamo nell’estensione PySparse, in tre formati differenti. Il primo, LL, e piu versatile nel

creare e manipolare la matrice sparsa, mentre gli altri due, Compressed Sparse Row (CSR)

e Sparse Symmetric Skyline (SSS ), sono piu efficienti nell’eseguire con essa operazioni di

tipo BLAS2 o BLAS3. E’ preferibile dunque convertire la matrice in un formato efficiente

prima di eseguire tali operazioni. A questo scopo sono disponibili le funzioni opportune.

NB: ricordiamo che in una matrice sparsa ci sono in genere molte componenti non im-

plementate, e trattate come zeri, e relativamente poche componenti implementate, dette

nonzeri. Le prime sono zeri strutturali, le seconde sono valori in generale diversi da zero

ma che possono anche essere uguali a zero (zeri accidentali). La disposizione dei nonzeri si

17


dice anche struttura di sparsita della matrice. Il rapporto tra numero di nonzeri e numero

totale di componenti ’e detto indice di densita, ed e compreso nell’intervallo [0, 1]. Una

matrice con indice di densita pari a 1 si dice densa.

Vediamo una sintesi delle operazioni principali:

from spmatrix import * # caricamento del modulo ’spmatrix’

# dell’estensione PySparse

from spmatrix_util import * # caricamento del modulo ’spmatrix_util’

# dell’estensione PySparse

• Creazione di Matrici Sparse (in formato LL):

M = ll_mat(nrighe, ncol) # ogni componente e’ di tipo ’double’

M = ll_mat(nrighe, ncol, nnz) # cosi’ viene inizialmente allocato spazio

# per ’nnz’ componenti (di default sono 1000)

M.shape # restituisce una tupla con la lunghezza di ogni dimensione di "v"

M.nnz # questo attributo contiene il numero di ’nonzeri’:

0 # appena creata, una matrice sparsa non ha ’nonzeri’

float(M.nnz)/(M.shape[0]*M.shape[1]) # indice di densita’ di M

M = ll_mat_rand(nrighe, ncol, ids) # genera una matrice sparsa e la inizializza

# con valori random e posizione casuale dei

# ’nonzeri’ scelta in modo che l’indice di

# densita’ risultante sia minore o uguale

# a ’ids’. I valori stanno in [0.0, 1.0)

B = M.copy() # crea una copia di ’M’

• Conversione di formato:

M.to_csr() # conversione di ’M’ al formato CSR

M.to_sss() # conversione di ’M’ al formato SSS

• Indicizzazione di Vettori e Matrici:

l’accesso e la modifica di singole componenti della matrice sono disponibili in tutti

e tre i formati, ed hanno modalita analoghe a quelle viste per gli array di Numeric

(PySparse e costruita ”sopra” Numeric):

print M[0,1] # stampa il valore dell’elemento di indici 0,1

M[r,c] = 9.77 # assegna l’elemento di indici r,c al valore 9.77

4.2.1 Operazioni implementate in modo dedicato per il formato LL


lunghezza n):

b = matvec(A, v) # equivale a "b = A * v"

b = matvec_transp(A, v) # equivale a "b = transpose(A) * v"

18



b[i] = A[i,0]*v[0] + A[i,1]*v[1] + ... + A[i,n-1]*v[n-1]

• copia:

W = M.copy() # crea una copia di "M" e "W" ne e’ il riferimento

• calcolo della norma:

M.norm(p) # "p" \’e una stringa; valori possibili: ’1’,’2’,...,’inf’,’fro’

4.2.2 Operazioni implementate in modo ottimizzato per il formato CSR ed SSS


lunghezza n):

b = matvec(A, v) # equivale a "b = A * v"

b = matvec_transp(A, v) # equivale a "b = transpose(A) * v"


b[i] = A[i,0]*v[0] + A[i,1]*v[1] + ... + A[i,n-1]*v[n-1]

5 Creazione di Grafici

Utilizziamo il modulo ’gplt’ facente parte dell’estensione SciPy. Esso si appoggia a Gnu-

plot.

import scipy # caricamento del modulo ’scipy’

from scipy import gplt # caricamento del modulo ’gplt’


5.1 Grafici bi-dimensionali

# caso generale: gplt.plot(x1,y1,format_string1, x2,y2,format_string2, ...)

# in cui tutti gli argomenti tranne ’y1’ sono opzionali

# e dove:

# ’xi’ e’ un array o una lista di ascisse (i=1,2,...)

# ’yi’ e’ un array o una lista di ordinate (i=1,2,...)

# ’format_string1,...2,...’ sono stringhe di formattazione

# che seguono lo stile di Gnuplot:

# ’title \’<title>\’ with <style>’, dove

# <title> puo’ essere un mnemonico che per il vettore,

# ad esempio ’dati misurati’, oppure ’notitle’;

# <style> puo’ essere ‘lines‘ (linea continua), ‘points‘

# (punti singoli), ‘linespoints‘, ‘impulses‘,

19


# ‘dots‘, ‘steps‘, ‘fsteps‘, ‘histeps‘, ...

# es. ’notitle with lines’

# esempi:

gplt.plot((1,2,3)) # grafica una lista di valori in funzione dell’indice,

# interpolandoli linearmente

gplt.plot([1,5,3]) # come sopra, con una lista

gplt.plot(array([1,2,3,5])) # come sopra, con un array

# caso completo:

x = arange(0.0, 1.0, 0.05) # scelgo un insieme di valori per le ascisse

y = x**3 # scelgo un insieme di valori per le ordinate

gplt.plot(x,y) # stampa per punti linearmente interpolati il grafico

# di ’y’ in funzione di ’x’

gplt.title(’grafico di funzione cubica’)

gplt.xtitle(’ascisse’) # stampa l’etichetta ’ascisse’ sull’asse delle ascisse

gplt.ytitle(’ordinate’) # stampa l’etichetta ’ordinate’ sull’asse delle ordinate

gplt.autoscale() # fa lo scalamento automatico del grafico

gplt.text((x,y),string) # stampa la stringa ’string’ alle coordinate ’x’,’y’

gplt.legend() # crea la legenda del grafico

gplt.hold(’on’) # i prossimi comandi ’plot’ vengono sovrapposti al

# grafico attuale

gplt.hold(’off’) # il prossimo comando ’plot’ sovrascrive il

# grafico attuale

Nota: quando si usa ”gplt.hold(’on’)” ricordarsi poi di mettere ”gplt.hold(’off’)” alla

fine, altrimenti ad ogni esecuzione successiva si accumulano i grafici delle precedenti

esecuzioni, inutilmente.

gplt.xaxis((min,max)) # visualizza l’intervallo dell’asse ’x’ tra ’min’ e ’max’

gplt.yaxis((min,max)) # visualizza l’intervallo dell’asse ’y’ tra ’min’ e ’max’

gplt.grid(’on’) # visualizza una griglia cartesiana

gplt.grid(’off’) # nasconde la griglia

gplt.logx(’on’) # imposta la scala logaritmica sull’asse ’x’

gplt.logx(’off’) # toglie la scala logaritmica sull’asse ’x’

# (ripristina la scala naturale)

gplt.logy(’on’) # imposta la scala logaritmica sull’asse ’y’

gplt.logy(’off’) # toglie la scala logaritmica sull’asse ’y’

# (ripristina la scala naturale)

gplt.close() # chiude il grafico corrente e distrugge la finestra corrispondente

5.2 Grafici tri-dimensionali

Per fare il grafico di una superficie, in 3D, definita per punti interpolati linearmente, si

usa la funzione

# Per fare il grafico di una superficie in 3D definita da un insieme di

20


# punti interpolati linearmente:

gplt.surf(x,y,z)

# in cui tutti gli argomenti tranne ’z’ sono opzionali

# e dove:

# ’x’ e’ un array o una lista di ascisse (i=1,2,...)

# ’y’ e’ un array o una lista di ordinate (i=1,2,...)

# ’z’ e’ una matrice di valori il cui elemento generico

# z[i,j] e’ relativo alle coordinate x[i] e y[j]

gplt.zaxis((min,max)) # visualizza l’intervallo dell’asse ’z’ tra ’min’ e ’max’

gplt.hidden(’remove’)

gplt.view(x_rotation, z_rotation, scale, z_scale) # uso:

# gli angoli sono espressi in gradi .

gplt.view("default") # "default"’ corrisponde a:

# x_rotation=60, z_rotation=30, scale=1, z_scale=1

gplt.text(x,y,z,string) # stampa la stringa ’string’ alle coordinate ’x’,’y’,’z’

# esempio:

x=(arange(50.)-25)/2.

y=(arange(50.)-25)/2.

r = sqrt(x[:,NewAxis]**2+y**2)

z = scipy.special.j0(r)

gplt.surf(x,y,z)

#

# per fare il grafico di un insieme di punti in 3-D:

#

from RandomArray import normal

points = normal(0,1.,(500,3))

gplt.plot3d(points,’notitle with points’)

points = normal(2,1.,(500,3))

gplt.hold(’on’)

gplt.plot3d(points,’notitle with points’)

gplt.ticlevel(0)

5.3 Immagini

21


Parte II Librerie

6 Libreria di Algebra Lineare Numerica

6.1 Interfaccia di Numeric alla LAPACK

from Numeric import * # caricamento del modulo principale dell’estensione Numeric

from LinearAlgebra import * # caricamento del modulo LinearAlgebra di Numeric

• soluzione di un sistema di equazioni lineari rappresentato mediante una matrice A

di dimensione n × n e vettore colonna v di lunghezza n:

x = solve_linear_equations(A, b) # "x" e’ un vettore colonna di lunghezza "n"

• soluzione di un sistema di equazioni lineari con un membro destro multiplo, ovvero

di ”p” sistemi di equazioni lineari aventi la stessa matrice di coefficienti e differenti

vettori di coefficienti a membro destro, rappresentati mediante una matrice A di

dimensioni n × n ed una matrice B di dimensioni n × p:

X = solve_linear_equations(A, B) # "X" e’ una una matrice di dimensioni "n x p"

6.2 Soluzione di sistemi lineari con matrice sparsa

L’estensione PySparse (vedi documentazione) contiene moduli che implementano metodi

diretti e metodi iterativi:

• superlu contiene funzioni per il pre-ordinamento di matrici sparse, per la fattoriz-

zazione LU e per la conseguente soluzione del sistema lineare;

• itsolvers contiene funzioni per la soluzione di sistemi lineari mediante metodi

iterativi, ad es. l’algoritmo del Gradiente Coniugato Precondizionato.

7 Libreria per la soluzione di Equazioni Nonlineari

8 Libreria per l’Interpolazione ed Approssimazione di funzioni

e di dati

import scipy # caricamento del modulo ’scipy’


22


8.1 Interpolazione polinomiale

Utilizziamo il modulo ’interpolate’ facente parte dell’estensione SciPy.

from scipy import interpolate # carica il modulo ’interpolate’

# definizione di un polinomio:

p = poly1d([1,3,4]) # definisce il polinomio: x^2 + 3*x + 4

# derivata di un polinomio:

d = p.deriv() # "d" e’ il polinomio uguale alla derivata di "p"

# valutazione di un polinomio:

x = [1.2,3.55,7.0,14.9] # insieme di valori

y = p(x) # "y" e’ il vettore di valori che assume "p" nei punti "x"

# interpolazione lineare a tratti di una funzione uni-variata:

x = [1,2,3,4,5]

y = [2.1,3.0,6.1,2.5,6.99]

f = interpolate.interp1d(x,y)

xnew=[1,1.5,2,2.5]

gplt.plot(xnew,f(xnew))

# interpolazione di una funzione uni-variata mediante ’splines’ cubiche:

x = [1,2,3,4,5]

y = [2.1,3.0,6.1,2.5,6.99]

s1 = interpolate.splrep(x,y,s=0.0)

xnew=[1,1.5,2,2.5]

ys1 = interpolate.splev(xnew,s1,der=0)

gplt.plot(xnew,ys1)

8.2 Approssimazione ai Minimi Quadrati Discreti

Utilizziamo il modulo ’optimize’ facente parte dell’estensione SciPy.

from scipy import optimize # carica il modulo ’optimize’

L’approssimazione ai minimi quadrati, cioe la soluzione la soluzione x al problema

min ‖f(x)‖22, si ottiene chiamando la funzione leastsq(...) (vedere istruzioni con

help(optimize.leastsq)).

La funzione f(x) puo anche essere lineare, ad esempio f(x) = A · x− b, con matrice A di

dimensione m × n e vettori colonna x di lunghezza n e b di lunghezza m.

23


9 Libreria per l’Integrazione Numerica

10 Libreria per la soluzione di Problemi ai Valori Iniziali

11 Libreria per la soluzione di Problemi ai Limiti

12 Libreria per la Geometria Computazionale

13 Programmazione ad Oggetti in linguaggio Python

13.1 Classi e Oggetti

Implementazione della classe:

• definizione di una classe: class < nome classe > : . Nelle righe successive, vengono

dichiarate le variabili ed i metodi della classe, indentandole.

• definizione di una classe derivata: class < nome classe > ( < classe padre > ):

. Nel caso in cui un metodo non sia presente, viene cercato ricorsivamente nella

gerarchia delle classi padre; questo avviene anche per il costruttore.

• definizione del costruttore: def init (self , < argomenti > ) :

• definizione del distruttore: def del (self , < argomenti > ) : . Lo si usa cosı:

del < nomeistanza > .

• definizione di un metodo: def < nome metodo >(self , < argomenti > ) :

• definizione di una variabile: self.< nome variabile > = < valore >

Nota: self e il riferimento all’istanza (oggetto) che sta chiamando il metodo della classe.

Nota: i metodi di una classe si comportano come una normale funzione Python; e possibile

assegnare valori di default nel seguente modo: < argomento >=< valore di default >.

Utilizzo della classe:

• creazione di un’istanza della classe: < variabile > = < nome classe >()

• accesso ai suoi metodi : < nome istanza >.< nome metodo >( < argomenti > )

oppure < nome classe >.< nome metodo >( < nome istanza > , < argomenti >

)

• accesso alle sue variabili : < nome istanza >.< nome variabile >

Ispezione di classi e oggetti:

• la funzione (built-in) type( < nome > ) ci dice se nome e ’class’ oppure ’instance’.

24


• la funzione (built-in) dir( < nome > ) ritorna una lista dei nomi delle proprieta e

dei metodi, se nome e una classe, ed una lista dei nomi delle variabili, se nome e

un’istanza.

• la proprieta bases di una classe, contiene la lista dei nomi delle classi da cui

questa e derivata.

• la proprieta dict di un’istanza, contiene la lista delle coppie < nome >:<

valore > relative alle sue variabili interne.

• la proprieta class di un’istanza, contiene il nome della classe che l’ha generata.

13.1.1 Esempio: la classe Matrice

25

Calcolo Numerico in Python

Documents

Transcript of Calcolo Numerico in Python