Introduzione_Python

7/23/2019 Introduzione_Python

1/25

Calcolo Numerico in Python

Fabio Marcuzzi, Maria Rosaria Russo

Dipartimento di Matematica Pura ed Applicata

Universita di PadovaItalyemail: [email protected], [email protected]

Versione 3.0

6 Giugno 2005


2/25

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 decisioneif-then-else . . . . . . . . . . . . . . . . . . . . . . . . 10

3.2.2 Il ciclofor

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 144.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


3/25


5 Creazione di Grafici 19

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

5.2 Grafici tri-dimensionali . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205.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 lInterpolazione ed Approssimazione di funzioni e di dati 22

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

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

9 Libreria per lIntegrazione 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


4/25


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 dellutilizzo;

2. effettua in modo automatico lallocazione 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

Lesecuzione di programmi scritti in Python avviene tramite la sua shell interattiva, che

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

In generale, e utile dotarsi di unIDE (Integrated Development Environment, ovvero

un ambiente di sviluppo integrato). A questo proposito ci sono varie soluzioni. UnIDEdisponibile 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 lIDE 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


5/25


1.3 Moduli aggiuntivi per il calcolo scientifico

Lutilizzo di Python per limplementazione 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 lestensione 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

unestensione che offre vari strumenti di calcolo scientifico; si appoggia a Numeric.

PySparse ( ci riferiamo alla versione disponibile al sito

http://pysparse.sourceforge.net/ ): e lestensione che offre limplementazione

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 lestensione

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 unestensione alternativa a Pygist, che interfaccia Python

al popolare Gnuplot.

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

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

alla simulazione mediante il Metodo degli Elementi Finiti).

5


6/25


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 .py

2. invocare lambiente basilare di sviluppo IDLE e poi lanciare il programma:

idle

(viene aperta la finestra di IDLE)

selezionare File>Opene poi selezionare il file contenente il programma;

(viene aperta la finestra delleditor con il file selezionato)

sulla finestra delleditor con il programma da avviare, selezionare Run

>RunModule

(sulla shell di IDLE appaiono i risultati dellesecuzione e rimangono disponi-

bili).

2.1 Caricamento di Moduli aggiuntivi (import e reload)

In entrambe le modalita di avvio di un programma, precedentemente citate, lambiente di

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

disponibili, ma non sono stati inizializzati, con leccezione disys (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 listruzione import) ilmodulo 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 ()del generico modulo

, e possibile agire in uno dei seguenti modi:

caricamento di un modulo mantenendone il nome:

import

chiamata: .()

6


7/25


caricamento di un modulo nascondendone il nome:

from import *

chiamata: ()

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

from import [, , ...]

chiamata: ()

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 caricarlonuovamente nellambiente. In questo caso, pero, listruzione importnon farebbe nulla (in

quanto il modulo e gia stato importato), mentre e necessario utilizzare listruzionereload:

reload() .

2.2 Lavorare con la shell interattiva di Python

Dalla shell interattiva di Python e possibile impartire qualunque istruzione eseguibile

dallinterprete Python. Quindi: le istruzioni dei programmi possono essere singolarmente

provate sulla shell (ad es. per verificarne lesecuzione). 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 unistruzione precedentemente impartita alla shell di Python, nella

sessione corrente, puo essere facilmente richiamata:

con IDLE (Linux): muovendosi con le combinazioni di tasti e nellelenco storico progressivo memorizzato dallambiente;

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

combinazioni di tasti e nellelenco storico

progressivo memorizzato dallambiente .

Questa funzione e analoga a quella ottenuta in Matlab cone )

.

2.2.2 Interazione con le directories

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

import os

7


8/25


Nome della directory corrente:

os.getcwd()

Cambio della directory corrente:

os.chdir() # passa alla directory figlia specificata

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

os.chdir() # passa alla directory specificata

Elenco del contenuto di una directory:

os.listdir(.) # mostra la directory corrente

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

os.listdir() # mostra la directory figlia specificataos.listdir() # mostra la directory specificata

Creazione di una directory:

os.mkdir() # crea la directory figlia specificata

os.mkdir() # crea la directory specificata

2.2.3 Help in linea

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

un modulo intero:

help()

una funzione di un modulo:

help(.)

Nota: la funzione help() richiede che il 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


9/25


>>> 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


10/25


listino[2] = 3 # imposta a 3 lelemento di indice 2

del listino[2] # cancella lelemento 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 lindice dove viene trovata la

# prima occorrenza di "acqua"

list.append("vino") # aggiunge lelemento "vino" alla fine della lista

3.2 Istruzioni composte

3.2.1 La decisione if-then-else

Permette di esprimere unalternativa in base ad una condizione da testare durante

lesecuzione del programma :

if :

else:

#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 allinterno di un ciclo.

3.2.2 Il ciclo for

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

:

for in :

#end

10


11/25


Il puo essere espresso mediante una listaoppure un array (Nu-

meric.array), eassume 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


12/25


Lidentificativodef(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 lindentazione torna al livello della

funzione o quando incontra listruzione 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 allinterno della

funzione stessa, e di nascondere alle istruzioni di questa le altre variabili definite allesterno

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

#end

in_file.close() # chiusura del file

Nota: lestensione del file puo essere qualunque, e non necessariamente txt.

12


13/25


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 timevenga prima importato:

import time

dopodiche listante attuale, misurato in secondi, viene preferibilmente fornito dalla

funzione:

in Linux:

istante_1 = time.time() # listante attuale viene memorizzato

# nella variabile "istante_1"

in Windows:

istante_1 = time.clock() # listante attuale viene memorizzato

# nella variabile "istante_1"

Quindi, e molto semplice misurare lintervallo di tempo necessario allesecuzione 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


14/25


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

dovuto ad operazioni esterne allalgoritmo di calcolo numerico in esecuzione, e

che quindi prendere il tempo di esecuzione come tempo di calcolo pu o risultare

unapprossimazion 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

Lelaborazione di vettori e matrici non puo essere fatta con gli strumenti base di Python,

ad es.listee ciclifor, in quanto risulterebbe assai lenta se confrontata con implementazioniin linguaggio C o Fortran. Per questo motivo e stata sviluppata lestensioneNumeric, detta

anche NumPy, che consente di eseguire operazioni su vettori e matrici con unefficienza

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

La documentazione HTML relativa allestensione Numericpuo 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 doublecorrisponde al numero in virgola mobile in doppia precisione

(ildoubledel 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


15/25


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 2x3v = 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 lultimo 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 dellelemento di indici 0,1

M[r,c] = 9.77 # assegna lelemento 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 lesecuzione. 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 vv /= 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


16/25


Nota: lestensione 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 brisulta essere un vettore, di lunghezza n, il cui i-esimo componente e pari a:

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

prodotto per componentidei vettori v e w, di lunghezza n:

b = v * w


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

16


17/25



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 no 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 nellestensionePySparse, 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 nelleseguire 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 comezeri, 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 deinonzerisi

17


18/25


dice anchestruttura di sparsitadella matrice. Il rapporto tra numero di nonzerie numero

totale di componenti e detto indice di densita, ed e compreso nellintervallo [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

# dellestensione PySparse

from spmatrix_util import * # caricamento del modulo spmatrix_util

# dellestensione 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 lindice 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:

laccesso 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 dellelemento di indici 0,1

M[r,c] = 9.77 # assegna lelemento di indici r,c al valore 9.77

4.2.1 Operazioni implementate in modo dedicato per il formato LL

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

lunghezza n):

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

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

18


19/25



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

prodotto matrice-vettore (matrice A di dimensione n n e vettore colonna v dilunghezza 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 dellestensione 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 \\ with , dove

# puo essere un mnemonico che per il vettore,

# ad esempio dati misurati, oppure notitle;

# puo essere lines (linea continua), points

# (punti singoli), linespoints, impulses,

19


20/25


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

# es. notitle with lines

# esempi:

gplt.plot((1,2,3)) # grafica una lista di valori in funzione dellindice,# 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 letichetta ascisse sullasse delle ascisse

gplt.ytitle(ordinate) # stampa letichetta ordinate sullasse 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 lintervallo dellasse x tra min e max

gplt.yaxis((min,max)) # visualizza lintervallo dellasse 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 sullasse x

gplt.logx(off) # toglie la scala logaritmica sullasse x

# (ripristina la scala naturale)

gplt.logy(on) # imposta la scala logaritmica sullasse ygplt.logy(off) # toglie la scala logaritmica sullasse 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


21/25


# 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 lintervallo dellasse 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


22/25


Parte II Librerie

6 Libreria di Algebra Lineare Numerica

6.1 Interfaccia di Numeric alla LAPACK

from Numeric import * # caricamento del modulo principale dellestensione Nu

from LinearAlgebra import * # caricamento del modulo LinearAlgebra di Numeri

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

di dimensionen ne 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 ned 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

Lestensione 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. lalgoritmo del Gradiente Coniugato Precondizionato.

7 Libreria per la soluzione di Equazioni Nonlineari

8 Libreria per lInterpolazione ed Approssimazione di funzioni

e di dati

import scipy # caricamento del modulo scipy


22


23/25


8.1 Interpolazione polinomiale

Utilizziamo il modulo interpolate facente parte dellestensione 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 dellestensione SciPy.

from scipy import optimize # carica il modulo optimize

Lapprossimazione 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) =Axb, con matrice Adi

dimensione m ne vettori colonna xdi lunghezza n e bdi lunghezza m.

23


24/25


9 Libreria per lIntegrazione 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 allistanza (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 unistanza 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 nomeeclass oppure instance.

24


25/25


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

unistanza.

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

questa e derivata.

la proprieta dict di unistanza, contiene la lista delle coppie < nome >: relative alle sue variabili interne.

la proprieta class di unistanza, contiene il nome della classe che lha generata.

13.1.1 Esempio: la classe Matrice

Introduzione_Python

Documents

Transcript of Introduzione_Python