Appunti Programmazione - calcolo numerico - Python

Appunti di programmazione

Esame di calcolo numerico

Docente: Francesca Pitolli

Università degli Studi di Roma La Sapienza

Ingegneria Meccanica

Programmazione con Python

Costrutto if (if then else) - esempi

Esegue in alternativa diversi blocchi di istruzioni, una alternativa esclude l'altra. Se non è verificato il primo blocco passa al secondo, se non è verificato neanche il secondo passa al terzo e così via. Se è verificato il primo non passa a verificare il secondo! Attenzione all'ordine dei blocchi!

c=4
# c%4 mi dà il resto della divisione fra c e 4
if c%3!=0:
    print "c non è divisibile per 3!"
elif c%3==0:
    print "c è divisibile per 3!"

b=12
if b%3!=0:
    print "b non è divisibile per 3!"
else:
    print "b è divisibile per 3!"

d=7
if d%3==0 and d%2==0:
    print "d è divisibile per 3 e per 2"
elif d%3==0:
    print "d è divisibile per 3"
elif d%2==0:
    print "d è divisibile per 2"
else:
    print "d non è divisibile né per 2 né per 3"

Costrutto try-except - esempi

È un sistema di controllo dell'errore! È molto presente nelle librerie come numpy e math.

import sys  # è una libreria che lavora sul sistema operativo!

try:
    val=raw_input("Inserisci un numero")
    print float(val), '\n'
except:
    print "Devi inserire un numero!"
    sys.exit()

f=lambda t: -1.0/t  # è una funzione che voglio definire. la lambda function è un altro modo per definire una funzione diverso da def

try:
    X=raw_input("Inserisci un numero")
    X= float(X)
    print f(X)
except:
    print "Valore impossibile!"

Spiegazione: try indica i comandi che voglio fare, except indica cosa deve fare il programma se ciò che ho chiesto nel try non è corretto. Nel secondo esempio col try chiedo al programma di inserire dei dati di input (raw_input) e poi trasformare quel dato in float e poi stampare la funzione f (definita prima come lambda function) calcolata in quel valore. Se inserisco un valore che non ha senso per la funzione, come 0, allora l'except fa sì che il programma mi dica che il valore è impossibile e si blocchi.

Costrutto while - esempi

i=0 # condizione iniziale con cui entro nel ciclo while (deve soddisfare la condizione)

while i<10:  # condizione fino alla quale voglio che sia ripetuto il mio ciclo di istruzioni
    i=i+1  # ciclo di istruzioni che voglio siano ripetute fino a che non è verificata la condizione imposta
    print i

while 0<i<50:
    i=i+1
    print i

b=1000000
while b>=10:
    print b, b>=0
    b=b/10

A volte il costrutto while si può sostituire al for ma ha un costo computazionale più alto. Non conviene usarlo in questi casi. Dice al programma di ripetere un certo ciclo di istruzioni finché non è verificata una certa condizione. Si usa per i metodi iterativi.

Il modulo numpy - appunti

import numpy as np  # numpy è una libreria di calcolo numerico

# dir(numpy) per avere le funzioni di numpy e help(numpy.) per avere informazioni sulla singola funzione.
# per esempio help(numpy.random) mi dà le info su cosa fa numpy.random.
# oppure numpy.linalg (calcola il determinante di una matrice, il rango e altre cose simili sulle matrici.
# numpy consente di creare array numerici, vettori numerici o matrici (array a più indici)

# arange --> come la xrange
B=np.arange(10)  # vettore che va da 1 a 9
# posso trasformare una lista in un array di numpy con np.array(C)
C=[i for i in range(10)]
np.array(C)
# gli elementi del vettore si richiamano esattamente come nelle liste. esempio:
B[1:-1]
# se faccio dir(B) ho diversi attributi che per le liste
B*2  # è diverso moltiplicare una lista o un vettore per un numero. B*2 è diverso da C*2 perché B è un vettore e C una lista. C viene ripetuta come lista mentre B ogni elemento viene moltiplicato per il numero. idem per le altre operazioni algebriche

D=B.reshape()  # cambia ordine degli elementi ottenendo per esempio una matrice.
# esempio: D=B.reshape(5,2)
D[1,1]  # estraggo l'elemento della seconda riga e seconda colonna
D[0,0]  # estraggo elemento prima riga e prima colonna
D[:,1]  # estraggo tutte le righe della seconda colonna
D[:,1]=np.array([11,12,13,14,15])  # cambio la seconda colonna di D. oss: abbiamo modificato anche B perché D era il reshape per B cioè ho fatto una copia per reference, cioè ho lavorato sulla stessa zona di memoria. non ho fatto una vera copia in un'altra zona di memoria (per cui dovrei usare B.copy())

B.sum()  # somma gli elementi del vettore
B.prod()  # moltiplica gli elementi di B
B.cumsum  # somma cumulativa degli elementi
D.sum(axis=0)  # sommo la prima colonna
D.sum(axis=1)  # sommo la riga
# con gli array posso fare le operazioni sul vettore senza dover estrarre elemento per elemento e fare l'operazione
# posso fare anche funzioni di vettori. non c'è bisogno di usare la libreria math perché numpy le ha già al suo interno

np.sin(D)  # seno di ogni elemento di D
np.max(D,axis=0)  # massimo su ogni riga
np.mean(D)  # media su tutti gli elementi
np.sort(np.abs(B))  # riordino gli elementi in valore assoluto (sort per riordinarli e abs per fare il valore assoluto)
# posso fare dei confronti: esempio: 2>3 mi dà come output falso, se faccio 2>=3 mi dà come output vero

B>np.mean(B)  # quali elementi di B sono maggiori della media degli elementi di B?
# con questi confronti posso creare una maschera che filtri, estragga alcuni elementi

V=np.random.random_integers(0,1,size=10)
np.nonzero(V)  # estrae gli elementi non nulli del vettore V. richiama solo gli elementi di V non nulli --> crea una maschera
A=np.reshape([i for i in range(1,10)],(3,3))  # faccio la lista da 1 a 9 e poi con reshape la riordino in forma 3x3
# se definisco una matrice A posso fare la trasposta con : A.T
# posso mettere insieme comandi e stringhe:
print "A:", A, "\n A trasposta:", A.T
np.diag(A)  # estraggo la diagonale principale
np.diag(A,-2)  # estraggo la diagonale -2
np.tril(A)  # estrae la parte triangolare superiore
np.triu(A)  # estrae la parte triangolare inferiore
np.triu(A,-1)  # estrae parte triangolare superiore a partire dalla diagonale -2
np.trace(A)  # traccia della matrice A

np.dot(A,L)  # prodotto righe per colonne fra due matrici A e L (se sono vettori prodotto scalare)
# attenzione alle dimensioni delle matrici: valgono le regole di algebra lineare.
np.dot(V,V) = np.sum(V**2)

A*V  # se faccio la moltiplicazione e basta? A è una matrice 3x3 e V un vettore con 3 elementi. trasforma V in una matrice che ha tutta la prima colonna col primo elemento, tutta la seconda col secondo elemento, tutta la terza col terzo e poi moltiplica. fa il cosiddetto broadcast.

Lettura dati di input come stringa - esempio

(Esempio coi dati di input per il metodo upwind)

import sys

# Lettura dati di input (come stringa)
dati_input = raw_input("Inserisci i valori c T0 T a b M N separati da uno spazio: ")

# La raw_input è una funzione che consente di inserire dei dati di input come stringa (da manipolare poi) richiedendoli direttamente con istruzioni

# Trasforma la stringa di input in numeri
dati = []  # crea un elenco di nome dati
for valore in dati_input.split():  # spacchetta i dati della raw_input in prossimità degli spazi
    try:
        dati.append(float(valore))  # aggiunge gli elementi alla lista dati come float
    except:  # comando di controllo
        print("I numeri devono essere separati da uno spazio!!")
        sys.exit()

if len(dati) != 7:  # comando di controllo
    print("Devi inserire 7 numeri!!")
    sys.exit()

# Dati di input del problema
c, T0, T, a, b, M, N = dati[0], dati[1], dati[2], dati[3], dati[4], int(dati[5]), int(dati[6])  # M e N devono essere interi perché stanno nei contatori
print "c: ", c, "T0: ", T0, "T: ", T, "a: ", a, "b: ", b
print "M: ", M, "N: ", N
# Con sistemi di controllo per evitare errori
# Meglio usare questo modo di immissione dati di input perché così definisco direttamente quali dati devono essere int e quali float ed evito errori

Disegnare e dare i valori della funzione coseno

import math
import matplotlib.pyplot as plt

nodi = [i*2.0*math.pi/20 for i in range(21)]
valori = [math.cos(nodi[i]) for i in range(21)]

print nodi
print valori

plt.figure()
plt.plot(nodi, valori)
plt.show()

Metodo di Eulero

Problema di Cauchy:

• y' = -cosy^2 + sint^2
• y(0) = sqrt(2)
• Risoluzione con metodo di Eulero

# Librerie importate:
import math  # libreria che contiene le funzioni matematiche elementari
import matplotlib.pyplot as plt  # libreria dei grafici. "as plt" significa che ho "soprannominato" questa libreria "plt" per semplicità

# Definizione della funzione termine noto f(t,y)
def f(t,y):
    return -math.cos(y**2)+math.sin(t**2)  # davanti alle funzioni elementari devo sempre mettere math. perché devo dire al programma che deve pescare quella funzione dalla libreria math

# Discretizzazione
beta = 4.0  # ampiezza intervallo discretizzazione
n = 8  # numero di nodi. se aumento il numero di nodi (8, 16, 32,40...) la funzione è approssimata sempre meglio.
h = beta/n  # definizione passo discretizz