Pel setembre d’aquest any passat vaig estar per Barcelona a la PHPConference’09, i allà hi havia un estant de l’editorial O’Relly. D’aquestes coses que un es possa a veure titols de llibres (la veritat tenen títols molts interessants) en vaig trobar un sobre concurrència que la veritat tenia molt bona pinta, de la serie Theory/In/Practice anomenat The Art of Concurrency.
Això que el mateix vespre en vaig possar a llegir-lo, i hem vaig trobar un exemple com calcular el valor de la constant 
emprant integració numèrica. Per a capitols posteriors mostrar de com descomposar el problema i paral·lelitzar-lo. El càlcul de la constant ho fa integrant numèricamet entre ![[0,1]](http://s.wordpress.com/latex.php?latex=%5B0%2C1%5D&bg=T&fg=000000&s=0 "[0,1]")
la funció 
, dins el meu infim coneixement sobre anàlisi i mètodes numèrics vaig mirar analìticament de fer la integral definida:
![\pi=\int ^1 _0 \frac{4}{1+x^2}\ dx=4 \int ^1 _0 \frac{dx}{1+x^2}=4[\arctan 1-\arctan 0]\thickapprox 3.1415926535897931](http://s.wordpress.com/latex.php?latex=%5Cpi%3D%5Cint%20%5E1%20_0%20%5Cfrac%7B4%7D%7B1%2Bx%5E2%7D%5C%20dx%3D4%20%5Cint%20%5E1%20_0%20%5Cfrac%7Bdx%7D%7B1%2Bx%5E2%7D%3D4%5B%5Carctan%201-%5Carctan%200%5D%5Cthickapprox%203.1415926535897931&bg=T&fg=000000&s=1 "\pi=\int ^1 _0 \frac{4}{1+x^2}\ dx=4 \int ^1 _0 \frac{dx}{1+x^2}=4[\arctan 1-\arctan 0]\thickapprox 3.1415926535897931")
Com sabem l’ordinador no resol a nivell analític si no que tot és numèric, funcions com la del 
i d’altres es fan mitjançant eines numèriques com les series de MacLaurin per fer aproximacions. En el nostre cas, seria obtenir el valor aproximat de emprant integració nunèrica, més concretament la regla del rectàngle que és una de les regles de quadratura numèrica que s’ensenyen durant la carrera. En aquell moment no m’hi vaig fixar, però, per què l’autor va triar aquest mètode com exemple? Tots sabem que el mètode en integració numèrica que convergeix més rapidament és el de Simpson i el que en poques iteracions obtenim el resultat esperat amb el mínim error. Si miram el codi d’exemple que possa ell, emprava unes 100000 iteracions per obtenir un valor molt aprop de
-
float calcPiRect(long num_rects)
-
{
-
int i;
-
double mid, height, width, area, sum=0.0;
-
-
width=1.0/(double) num_rects;
-
for (i=0;i<num_rects ;i++)
-
{
-
mid=(i+0.5)*width;
-
height=4.0/(1.0+mid*mid);
-
sum+=height;
-
}
-
return(sum*width);
-
}
Fa poc dies vaig tornar a agafar el llibre i vaig fer la versió però emprant Simpson, i ja que hi era amb Python. Vaig fer una petita variació, en aquesta versió. Empro una variable pels valors dels nodes 
en lloc d’emprar un array, això sacrificaria una mica de memòria. La veritat vaig notar la diferència només fent 14 iteracions s’aproximava molt al valor real amb un error molt baix (
) en front a l’error donat per l’algorisme de rectàngle amb les mateixes iteracions (
)
-
def func(x):
-
return(4.0/(1+x*x))
-
-
def simpson(a, b, n):
-
if n< =0:
-
raise Exception, "n ha d'esser mes gros que 0"
-
-
if a>=b:
-
raise Exception, "a ha d'esser mes petit que b"
-
-
h = (float(b)-float(a))/n
-
ini = h/6
-
sum = 0.0
-
-
xn = a
-
for i in range(1, n+1):
-
xn += h
-
sum += (func(xn-h)+4*func((xn+(xn-h))/2)+func(xn))
-
-
return(ini*sum)
Com es pot observar a la taula i a la gràfica (en escala logarítmica) següents el mètode del rectangle és molt lent a convergir cap a , en canvi el de Simpson amb 14 iteracions ja tenim un valor amb 10 decimals estable amb un error més que acceptable. Per això moltes vegades i en aquests tipus d’algorismes no només s’ha de veure el cost computacional si no que també la velocitat de convergència, ja que això ens reduirà el temps de càlcul del problema en el que estam treballant.
| # |
M. Rectangle |
M. Simpson |
 Rectangle |
Simpson |
| 1 |
0.0000000000 |
3.1333333333 |
3.1415926536 |
0.0082593203 |
| 2 |
1.8823529412 |
3.1415686275 |
1.2592397124 |
0.0000240261 |
| 3 |
2.3639639640 |
3.1415917809 |
0.7776286896 |
0.0000008727
|
| 4 |
2.5804288370 |
3.1415925025 |
0.5611638166 |
0.0000001511 |
| 5 |
2.7029369137 |
3.1415926139 |
0.4386557399 |
0.0000000397 |
| 6 |
2.7816432763 |
3.1415926403 |
0.3599493773 |
0.0000000133 |
| 7 |
2.8364440025 |
3.1415926483 |
0.3051486511 |
0.0000000053 |
| 8 |
2.8767824635 |
3.1415926512 |
0.2648101901 |
0.0000000024 |
| 9 |
2.9077111792 |
3.1415926524 |
0.2338814744 |
0.0000000012 |
| 10 |
2.9321763135 |
3.1415926530 |
0.2094163401 |
0.0000000006 |
| 11 |
2.9520110875 |
3.1415926532 |
0.1895815661 |
0.0000000003 |
| 12 |
2.9684157005 |
3.1415926534 |
0.1731769530 |
0.0000000002 |
| 13 |
2.9822087322 |
3.1415926535 |
0.1593839214 |
0.0000000001 |
| 14 |
2.9939675921 |
3.1415926535 |
0.1476250615 |
0.0000000001 |
| 15 |
3.0041112145 |
3.1415926535 |
0.1374814391 |
0.0000000001 |
| 16 |
3.0129509199 |
3.1415926536 |
0.1286417337 |
0.0000000000 |
| 17 |
3.0207228750 |
3.1415926536 |
0.1208697786 |
0.0000000000 |
Share on Facebook
Ho he estat provant i he trobat una cosa curiosa. Per molt que puig el nombre d’itreracions (milers o milions), la precisió no puja. De fet, va oscil.lant.
Deu ser per la precisió de les operación en punt flotant?
Te pos l’algoritme (escrit un poc diferent) en Ruby:
def simpson(a, b, n)
h = (b.to_f – a.to_f) / n
sum = (1 … n * 2).inject(0) { |t, x| t += ((x % 2) + 1) * yield((a + h * x)/2) }
h / 6 * (yield(a) + 2 * sum + yield(b))
end
pi = simpson(0, 1, 16) { |x| 4.0 / (1 + x * x) }
Si és cert, es que inclús es pot calcular el limit en el qual permet saber quina serà la precisió màxima de l’algorisme.
Crec que he confirmat les sospites. Usant aritmetica “sencera” (emprant el tipus rational de ruby), a mesura que augment el nombre d’iteracions la precisió augmenta. Cosa que amb punt flotant no passa.
Amb aquests mètodes es pot saber també la fita d’errors, per exemple en Simpson tendriem que la fita de l’error és $latex \frac{(b-a)^5}{2880}K_5$ on $latex K_5$ és la derivada 5 de la f nostre funció que calculam.