Calcolo per l’Astronomia a.a. 2014/2015 

Lauro Moscardini lauro.moscardini@unibo.it 

Tel. 051‐2095726 Sito web:  

apps.difa.unibo.it/files/people/lauro.moscardini/calcolo_14‐15 

Programma del corso  

•  Integrazione di funzioni • Metodi di interpolazione di daG 

•  Generazione di numeri casuali e metodi Monte Carlo 

•  Equazioni differenziali ordinarie •  Equazioni differenziali parziali 

Integrazione di funzioni 

Integrazione di funzioni: caso unidimensionale 

Il problema generale e’: 

€

I = f (x)dxa

b

∫

Grazie al teorema fondamentale del Calcolo Integrale, si cerca la funzione primiBva F(x) tale per cui F’(x)= f(x). Dopo di che si ha che:  

€

I = f (x)dxa

b

∫ = F(x) ab = F(b) − F(a)

Formule di Integrazione di  Newton‐Cotes 

Strategia: sosGtuire una funzione complicata con una sua approssimazione facile da integrale: 

€

I = f (x)dxa

b

∫ ≈ fn (x)dxa

b

∫

dove fn(x) e’ un polinomio di forma 

€

fn x( ) = a0 + a1x1 + a2x

2 + a3x3 + ...+ an−1x

n−1 + anxn

Ordine 1: Regola del trapezoide In questo caso il polinomio e’ al prim’ordine (reOa): 

€

I = f (x)dxa

b

∫ ≈ f1(x)dxa

b

∫

€

f1 x( ) = a0 + a1x

€

f1 x( ) = f (a) +f (b) − f (a)b − a

x − a( )

Equazione della reUa passante per i due punG [a, f(a)] e [b, f(b)]  

Ordine 1: Regola del trapezoide Integrando il polinomio si ha: 

€

I ≅ f (a) +f (b) − f (a)b − a

x − a( )

a

b

∫ dx

€

I ≅ b − a( ) f (a) + f (b)2

SBma dell’errore totale locale: 

€

Et = −112

′ ′ f ξ( ) b − a( )3 dove ξ punto tra a e b 

CommenB:  ‐  l’errore cresce con il cubo della lunghezza dell’intervallino ‐  il risultato e’ esaUo se l’integranda e’ una funzione costante o lineare (f”=0)  

Ordine 1: Regola del trapezoide Come fare a migliorare?  Ridurre la lunghezza dell’intervallino, facendo una applicazione mulBpla del 

metodo su n intervallini! 

Commento:  ‐  l’integrale e’ approssimato dal prodoUo tra la lunghezza dell’intervallino e una sGma dell’altezza media della funzione sullo stesso intervallino  

€

h =b − an

€

I = f (x)dxa

a+h

∫ + f (x)dxa+h

a+2h

∫ + ...+ f (x)dxa+(n−2)h

a+(n−1)h

∫ + f (x)dxa+(n−1)h

b

∫

€

I ≅ b − a( )f (a) + 2 f (a + ih) + f (b)

i=1

n−1

∑2n 

Ordine 1: Regola del trapezoide Come varia la sBma dell’errore? L’intervallo e’ ora lungo  

Commento:  ‐  dimezzando la lunghezza degli intervallini, l’errore si riduce di un faUore 4! e cosi’ via… AUenzione pero’: a un certo punto si diventa dominaG dall’errore di arrotondamento!  

  €

Ei = −112

′ ′ f ξ i( ) b − an

3

€

h =b − an

Su ogni singolo intervallino si ha 

€

Et = −1

12n3b − a( )3 ′ ′ f ξ i( ) ≅

i=1

n

∑ −b − a( )3

12n2f "

Ordine 2: Regola di Simpson 1/3 In questo caso il polinomio e’ al second’ordine (parabola): 

€

I = f (x)dxa

b

∫ ≈ f2(x)dxa

b

∫

€

f2 x( ) = a0 + a1x + a2x2

€

h =b − a2

€

I ≅ b − a( ) f (a) + 4 f (a + h) + f (b)6

€

I ≅ h3f (a) + 4 f (a + h) + f (b)[ ]

Per definire una parabola servono  3 punB, quindi 

Ordine 2: Regola di Simpson 1/3 

SBma errore totale: 

€

Et = −h5

90f (4 ) ξ( ) = −

b − a( )5

2880f (4 ) ξ( )

dove ξ punto tra a e b 

Commento:  ‐  Nonostante sia al second’ordine, l’errore cresce con la quinta potenza della lunghezza dell’intervallino!!!  

€

h =b − an

€

I ≅ b − a( )f (a) + 4 f (a + ih) + 2 f (a + ih) + f (b)

i= 2,4,6

n−2

∑i=1,3,5

n−1

∑3n

Applicazione mulBpla su n intervallini 

€

Et ≅ −b − a( )5

180n4f (4 )

SBma errore totale si riduce a: 

Ordine 3: Regola di Simpson 3/8 In questo caso il polinomio e’ al terz’ordine (cubica): 

€

I = f (x)dxa

b

∫ ≈ f3(x)dxa

b

∫

€

f3 x( ) = a0 + a1x + a2x2 + a3x

3

€

h =b − a3

€

I ≅ b − a( ) f (a) + 3 f (a + h) + 3(a + 2h) + f (b)8

€

I ≅ 3h8

f (a) + 3 f (a + h) + 3 f (a + 2h) + f (b)[ ]

Per definire una parabola servono 4 punB, quindi 

Ordine 3: Regola di Simpson 3/8 

SBma errore totale: 

€

Et = −3h5

80f (4 ) ξ( ) = −

b − a( )5

6480f (4 ) ξ( )

dove ξ punto tra a e b 

Commento:  ‐  Stesso ordine di Simpson 1/3, ma con coefficiente minore   solo leggermente piu’ accurata (coefficiente numerico inferiore).  

Ordini superiori 

SBma errore totale: 

€

Et = −8h7

945f (6) ξ( ) dove ξ punto tra 

a e b 

Ordine 4: Regola di Boole 

€

I ≅ b − a( ) 7 f (a) + 32 f (a + h) +12(a + 2h) + 32 f (a + 3h) + 7 f (b)90

Ordine 5 

€

I ≅ b − a( )19 f (a) + 75 f (a + h) + 50(a + 2h) + 50 f (a + 3h) + 75 f (x + 4h) +19 f (b)288

dove ξ punto tra a e b 

SBma errore totale: 

€

Et = −275h7

12096f (6) ξ( )

Estrapolazione di Richardson 

•  L’estrapolazione di Richardson rappresenta un metodo in cui si combinano tra di loro due diverse sGme numeriche dell’integrale per oUenerne una piu’ accurata, sfruUando la conoscenza dell’andamento dell’errore. 

•  La sua implementazione numerica piu’ efficiente e’ il metodo di Romberg, una tecnica ricorsiva ideale per generare una sGma dell’integrale con una tolleranza sull’errore pre‐definita. 

Estrapolazione di Richardson 

Supponiamo di avere applicato un dato metodo per sGmare l’integrale. Avremo: 

€

I = I(h) + E(h) I(h) e’ la sBma dell’integrale E(h) e’ il suo errore di troncamento 

€

Et ≅ −b − a( )3

12n2f " = −

b − a12

h2 f "

€

h =b − an

Considerando il trapezoide con spaziatura h: 

N.B. Il trapezoide ha un’accuratezza O(h2) 

Estrapolazione di Richardson Supponiamo ora di avere oUenuto due sGme dell’integrale con 

due spaziature diverse, h1 e h2. Avremo: 

€

I = I(h1) + E(h1) = I(h2) + E(h2)Considerando indipendente da h il valore medio di f”, si ha:    

€

E h1( )E h2( )

=h12

h22

€

I(h1) + E h2( ) h1h2

2

= I(h2) + E(h2)€

E h1( ) = E h2( ) h12

h22 

sosGtuendo 

Estrapolazione di Richardson Da cui possiamo ricavare: 

€

E(h2) ≅I h1( ) − I h2( )1− h1 /h2( )2

Pertanto possiamo migliorare la sGma precedente:    

€

I = I(h2) + E(h2) ≅ I h2( ) +I h1( ) − I h2( )1− h1 /h2( )2

In parBcolare se h2=h1/2 

€

I ≅ I h2( ) +I h2( ) − I h1( )22 −1

=43I h2( ) − 1

3I h1( )

Anche se basata sul metodo del trapezoide (applicato due volte!), la nuova sBma ha un’accuratezza O(h4)! 

Estrapolazione di Richardson Il metodo puo’ essere ulteriormente iterato. Se si hanno due 

diverse sGme con accuratezza O(h4) oUenute sempre con h2=h1/2, queste possono essere combinate per oUenere una nuova sGma che avra’ accuratezza O(h6): 

€

I ≅ 1615I h2( ) − 1

15I h1( )

Combinando due sGme con accuratezza O(h6) si puo’ oUenere una nuova sGma che avra’ accuratezza O(h8): 

€

I ≅ 6463I h2( ) − 1

63I h1( )

e cosi’ via… 

Algoritmo di integrazione di Romberg Corrisponde all’implementazione numerica dell’estrapolazione di Richardson, 

sempre con h2=h1/2. In generale i coefficienG necessari per combinare le diverse sGme possono essere scrif come: 

€

I j,k =4 k−1I j+1,k−1 − I j ,k−1

4 k−1 −1k rappresenta l’ordine di accuratezza: k=1 O(h2), k=2 O(h4), k=3 O(h6), … 

j rappresenta il livello di accuratezza: (j+1) e’ piu’ accurato di j 

O(h2)      O(h4)         O(h6)         O(h8)        O(h10) 

Algoritmo di integrazione di Romberg: pseudo‐codice I 

INPUT a, b (limiG di integrazione), INPUT εs (criterio di stop),INPUT maxit (numero massimo di iterazioni) n =1 CALL TRAPEZOIDE (n,a,b,integrale) I1,1= integrale εa=1.1 εs  i=0 DOWHILE (εa> εs) and (i<maxit)      i=i+1      n=2i 

     CALL TRAPEZOIDE(n,a,b,integrale) 

        I i+1,1= integrale       DOFOR k=2 to i+1              j=2+i+k 

       ENDDO 

ENDDO 

€

I j,k =4 k−1I j+1,k−1 − I j ,k−1

4 k−1 −1

€

εa =I1,i+1 − I1,iI1,i+1

⋅100

Algoritmo di integrazione di Romberg: pseudo‐codice II (implementata!) 

INPUT a, b (limiG di integrazione), INPUT εs (criterio di stop),INPUT maxit (numero massimo di iterazioni) 

εa=1.1 εs  i=0 DOWHILE (εa> εs) and (i<maxit)      i=i+1      n=2i‐1 

     CALL TRAPEZOIDE(n,a,b,integrale) 

        I i,1= integrale       DOFOR k=2 to i              j=1+i‐k 

       ENDDO 

ENDDO 

€

I j,k =4 k−1I j+1,k−1 − I j ,k−1

4 k−1 −1

€

εa =I1,i − I1,i−1I1,i

⋅100

QUADRATURA DI GAUSS •  Nelle formule di Newton‐Cotes, sono uGlizzaG valori della funzione 

in punB predeterminaB regolari.  

•  Nella Quadratura di Gauss, i punG di integrazione sono indeterminaB. 

•  Nell’esempio, l’approssimazione dell’integrale puo’ essere migliorata usando due punB interni scelB in modo ‘intelligente’. 

•  La Quadratura di Gauss offre un metodo per trovare quesG punG. 

QUADRATURA DI GAUSS: metodo dei coefficienB indeterminaB Per il trapezoide avevamo: 

€

I ≅ b − a( ) f (a) + f (b)2

=b − a2

f a( ) +b − aa

f b( )

Generalizzando: 

€

I ≅ c0 f x0( ) + c1 f x1( )dove x0 e x1 sono punG incogniG interni all’intervallo di integrazione e c0 e c1 sono opportuni coefficienG, pure da determinarsi  abbiamo 4 incognite da determinare! 

QUADRATURA DI GAUSS: Servono quaUro condizioni per trovare le quaUro incognite: 

€

c0 f x0( ) + c1 f x1( ) = 1dx−1

1

∫ = 2

c0 f x0( ) + c1 f x1( ) = xdx−1

1

∫ = 0

c0 f x0( ) + c1 f x1( ) = x 2dx−1

1

∫ =23

c0 f x0( ) + c1 f x1( ) = x 3dx−1

1

∫ = 0

NOTA: con un semplice cambiamento di variabile, qualsiasi integrale puo’ essere trasformato in un integrale tra ‐1 e 1! 

QUADRATURA DI GAUSS: 

SosGtuendo si ha: 

€

c0 + c1 = 2c0x0 + c1x1 = 0

c0x02 + c1x1

2 =23

c0x03 + c1x1

3 = 0

NOTA: con la scelta delle condizioni imposte, i valori risultanG delle incognite (c0, c1, x0, x1) sono tali da garanGre che l’integrale sia esaOo fino all’ordine cubico! 

QUADRATURA DI GAUSS: 

Risolvendo si ha: 

€

c0 =1c1 =1

x0 = −13≈ −0.5773503

x1 =13≈ 0.5773503

Pertanto: 

€

I ≅ f − 13

+ f 1

3

Formula a due punB di 

Gauss‐Legendre 

Per cambiare i limiB di integrazione, e’ necessario fare il seguente cambio di variabile: 

€

x = a0 + a1xdDall’estremo inferiore (x=a  xd=‐1) 

Dall’estremo superiore (x=b  xd=1) 

€

a = a0 + a1 ⋅ −1( )

€

b = a0 + a1 ⋅ +1( )

Le soluzioni sono: 

€

a0 =b + a2

€

a1 =b − a2

€

x =b + a( ) + b − a( )xd

2

€

dx =b − a( )2

dxd

Per cui 

il cui differenziale e’ 

QUADRATURA DI GAUSS: Formule agli ordini superiori 

In maniera analoga si possono sviluppare versioni della quadratura di Gauss agli ordini superiori. UGlizzando n punG si ha : 

€

I ≅ c0 f x0( ) + c1 f x1( ) + ...+ cn−1 f (xn−1)

Per esempio, la formula di Gauss‐Legendre a 3 punB ha i seguenG 6 parametri: 

€

c0 = 0.555555556c1 = 0.888888889c2 = 0.555555556

€

x0 = −0.774596669x1 = 0.0x2 = 0.774596669

Per oUenerli sono staG uGlizzaG 6 condizioni, pertanto la formula e’ esaOa per funzioni fino all’ordine 5 compreso!    

€

Et ≈ f (6) ξ( )

Formule agli ordini superiori 

€

Et ≈ f (6) ξ( )

€

Et ≈ f (4 ) ξ( )

€

Et ≈ f (8) ξ( )

€

Et ≈ f (10) ξ( )

INTEGRALI IMPROPRI Bisogna fare un cambio di variabile! 

€

f x( )a

b

∫ dx =1t 21/ b

1/ a

∫ f 1t

dt

Funziona quando f tende a zero piu’ velocemente di 1/x2 quando x tende a ∞ e per ab>0. AltrimenG conviene spezzare l’integrale in due parG. Esempio: 

€

f x( )−∞

b>0

∫ dx = f x( )−∞

−A

∫ dx + f x( )−A

b

∫ dx

FaUo normalmente FaUo cambiando variabile