Ecco un paio di modi per approssimare il valore di .
Partiamo da un classico, ma non troppo, che è detta serie di Leibniz-Gregory-Madhava (a seconda delle latitudini e delle epoche). Ovvero lo sviluppo in serie di potenze della funzione , per cui si ha che:
(1)
Che possiamo calcolare e stampare con il seguente codice OCTAVE.
printf('Approssimazione n° 1\n'); N = input('Sommare fino a N = '); piapp1 = 0; figure(1) for k=0:N piapp1 = piapp1 + 4*((-1)^k)/(2*k+1); subplot (211); hold on plot(k,piapp1,'*'); end plot(0:0.1:N,pi*ones(1,length(0:0.1:N)),'k--'); ylabel('Valore approssimato di \pi'); xlabel('Numero di termini della serie'); title(sprintf('Approssimazione n° 1, N = %d',N)); hold off
Che ha una convergenza molto, molto, lenta … per vederlo servono i numeri di Eulero, ma si può comunque fare qualche esperimento con il codice (un po’ goffo sotto molti punti di vista) riportato sopra.
Possiamo scegliere un’altra via, se sappiamo che
(2)
abbiamo che
(3)
Cosa ce ne facciamo? Be’ possiamo approssimare l’integrale con la regola dei trapezi … Ovvero discretizziamo l’intervallo in
punti equispaziati
, ovvero fissiamo una griglia su
di ampiezza
, per cui l’approssimazione dell’integrale diventa:
(4)
Si lo so, bisognerebbe pure sostituire quelle radici con un’approssimazione, ma facciamo finta di averlo fatto, o di averlo delegato ad OCTAVE.
printf('Approssimazione n° 2\n'); M = input('Numero dei punti in cui suddividere [0,1], M = '); for h=1:M x = 0:1/h:1; y = sqrt(1-x.^2); piapp2 = 4*trapz(x,y); subplot (212); hold on plot(h,piapp2,'*'); end plot(0:0.1:M,pi*ones(1,length(0:0.1:M)),'k--'); ylabel('Valore approssimato di \pi'); xlabel('Numero di intervalli di [0,1]'); title(sprintf('Approssimazione n° 2, M = 1,...,%d',M)); hold off
Il risultato dei due codici precedenti è questo:
Ah già, buon –day!
Codice dell’esempio: pigreco.m.tar.gz.
L’immagine di copertina viene dall’illustrazione all’inizio degli Elementi di Euclide nella traduzione di Adelardo di Bath (fonte: Wikipedia/Wikimedia).