Tag: pi greco

Buon 3/14 alle 9.26.53 (PM)!

Ecco un paio di modi per approssimare il valore di \pi.

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 \arctan(x), per cui si ha che:

(1)   \begin{equation*} \pi = 4 \sum_{k=0}^{+\infty} \frac{(-1)^k}{2k+1}. \end{equation*}

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)   \begin{equation*} \int \sqrt(1-t^2)dt = \frac{1}{2}\left(\sqrt{1-t^2}+\arcsin(t)\right) + C, \end{equation*}

abbiamo che

(3)   \begin{equation*} \pi = 4 \int_{0}^{1} \sqrt{1-t^2}dt. \end{equation*}

Cosa ce ne facciamo? Be’ possiamo approssimare l’integrale con la regola dei trapezi … Ovvero discretizziamo l’intervallo [0,1] in N+1 punti equispaziati \{x_k\}_{k=1}^{N+1}, ovvero fissiamo una griglia su [0,1] di ampiezza h = (1-0)/N, per cui l’approssimazione dell’integrale diventa:

(4)   \begin{equation*} \pi = 4 \int_{0}^{1} \sqrt{1-t^2}dt \approx \frac{1}{2N} \sum_{k=1}^N \left(\sqrt{1-(x_{k+1})^2}-\sqrt{1-(x_k)^2}\right). \end{equation*}

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 pi–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).