1 Trapezoidi e Simpson
Per prima cosa ci serve un array
che contenga i valori di una funzione in un certo intervallo, per esempio .
1 2 3 | N=100
x = linspace(0,2,N+1)
y = x**4 - 2*x + 1
|
Nella seconda riga sfruttiamo
la funzione linspace che genera un
array di valori equispaziati nell'intervallo . Nella seconda riga, sfruttiamo il fatto che ogni operazione
fatta su un array viene in automatico eseguita su ogni elemento di esso,
ottenendo cosi l'array della funzione che volevamo.
1 2 3 | h = abs( x[1] - x[0] )
trap = 0.5*h*(y[0] + y[-1])
trap += sum(y[1:-1]) * h
|
Definiamo la lunghezza
dell'intervallino di integrazione come differenza tra due valori consecutivi di
; scriviamo nella variabile trap la somma del primo e del'ultimo
valore della funzione divisi per due; successivamente, aggiungiamo al valore
appena registrato la somma di tutti i valori della funzione esclusi il primo e
l'ultimo.
Soffermiamoci per capire un pò
meglio alcune cose importanti, utili sempre, nascoste in queste apparentemente
innocenti righe:
- per avere un singolo elemento di un array usiamo x[i], con i la posizione dell'elemento partendo da zero; si possono usare anche indici negativi per partire dall'ultimo elemento dell'array, come abbiamo fatto in y[-1] per prendere appunto l'ultimo elemento di y;
- per poter scegliere alcuni elementi consecutivi di un array è possibile usare una slice : con x[i:j] creiamo un array con tutti gli elementi di x che vanno da quello di posto i a quello di posto j-1;
- mentre con trap = assegniamo alla variabile un valore (quindi se trap non esiste la creiamo, se contiene già un valore, questo viene sovrascritto), con trap += sommiamo al valore già contenuto in trap il valore che segue.
- la funzione abs(num) calcola il valore assoluto di num; la funzione sum(arr) restituisce la somma di tutti gli elementi del vettore arr
1 2 3 | simp = y[0] + y[-1]
simp += 4*sum(y[1:-1:2]) + 2*sum(y[2:-1:2])
simp *= h/3
|
- con y[1:-1:2] prendiamo dal secondo all'ultimo (escluso) saltando di due in due, quindi gli elementi di posto pari;
- con y[2:-1:2] prendiamo dal terzo all'ultimo (escluso) saltando di due in due, quindi gli elementi di posto dispari;
Torniamo al nostro problema:
l'integrale della funzione scelta. Dopo aver eseguito il codice visto sin ora,
dentro le variabili trap e simp abbiamo salvato il valore dell'integrale
calcolato con i due metodi. Poichè conosciamo il valore esatto dell'integrale
pari a , possiamo fare una prima verifica della bontà del nostro codice
e del metodo usato:
1 2 | esatto = 4.4
print "Valore esatto, Trapezoidi,
Simpson: ", esatto, trap, simp
|
Vediamo ora come fare un
analisi più estesa dell'andamento dell'errore dei due metodi in funzione del
numero di intervallini. L'idea di base è eseguire le righe di codice
precedenti per valori di compresi in un certo intervallo di valori. Ogni volta che un
nuovo valore dell'integrale viene calcolato, salviamo in un elemento di un array
l'errore, per entrambe i metodi. Questo lo si può fare in automatico utilizzando
un ciclo for, in questo modo:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | esatto = 4.4
Nmax=1000
trap_err=zeros((Nmax-1)/2.)
simp_err=zeros((Nmax-1)/2.)
i=0
for N in
range(2,Nmax,2):
x = linspace(0,2,N+1)
y = x**4 - 2*x + 1
h = abs( x[1] - x[0] )
trap = h* ( 0.5*y[0] + 0.5*y[-1] + sum(y[1:-1]) )
simp = h/3 *( y[0] + y[-1] + 4*sum(y[1:-1:2]) + 2*sum(y[2:-1:2]) )
trap_err[i] = abs(trap - esatto)
simp_err[i] = abs(simp - esatto)
i=i+1
|
- la funzione zeros(num) genera un array con num elementi di valore zero; trap_err e simp_err saranno dunque i vettori dove salveremo gli errori dei due metodi per ogni scelto;
- la funzione range(i,j,k) è spesso usata nei cicli for per generare una sequenza di valori nell'intervallo con passo
- il ciclo for esegue ripetutamente tutte le righe indentate e ad ogni ripetizione la variabile N assume i valori della sequenza fornita da range;
- il numero di volte che vengono eseguite è dato dalla lunghezza del sequenza fornita da range, pari a (Nmax-1)/2., valore con cui abbiamo dimensionato gli array che conterranno gli errori;
- in generale, per scrivere un valore in uno specifico elemento di posto di un array facciamo come nel codice sopra: trap_err[i] =.
Ora che abbiamo degli array
contenenti gli errori dei due metodi per un numero crescente di intervallini,
potremo farne il grafico in questo modo:
1 2 3 | N_arr = range(2,Nmax,2)
loglog(N_arr, trap_err)
loglog(N_arr, simp_err)
|
Creiamo prima un array con la
sequenza di valori usata nel ciclo for e poi usiamo la funzione loglog(x_arr,y_arr) che crea un grafico in
scala bi-logaritmica, con x_arr e
y_arr i vettori delle coordinate x,y
dei punti.
Il grafico che otteniamo è
del tutto analogo a quello mostrato a pag. (spiegare perché manca il rumore).
Come notiamo la pendenza delle due curve è diversa e dovrebbe darci l'ordine
dell'errore commesso con i due metodi. È quindi interessante verificarlo,
calcolandola in questo modo:
1 2 3 4 5 6 7 8 9 | deltaY = log(trap_err[350])-log(trap_err[300])
deltaX = log(N_arr[350])-log(N_arr[300])
slope_trap = deltaY / deltaX
deltaY = log(simp_err[350])-log(simp_err[300])
deltaX = log(N_arr[350])-log(N_arr[300])
slope_simp = deltaY / deltaX
print slope_simp, slope_trap
|
La pendenza la calcoliamo
geometricamente scegliendo opportunamente due punti per ogni grafico e facendo
il rapporto tra la differenza delle coordinate y e quella delle coordinate x.
Per esempio abbiamo scelto gli elementi 300 e 350 degli array trap_err, simp_err e N_arr e, poiché la relazione è lineare
in scala Log-Log, ne abbiamo fatto il logaritmo mediante la funzione log().