Sumar a cegues

L'entrada d'avui està dedicada a un petit incident de bar. És al bar on apareixen els problemes interessants.

Divendres, després de la V, em trobava al bar amb un company i vam demanar dos cafès amb llet i un croissant. A l'hora de pagar vaig demanar a última hora de pagar-ho tot junt. El resultat va ser que l'home del bar, em va demanar a mi de sumar-ho. Aquest va ser el meu procés mental:

1.04 + 1.04 + 1.24 = Uhm... uhh... 3x1 = 3 Llavors hi ha alguns cèntims que queden sueltos... Uhm... com a mínim 25 cèntims...

Mentre em trobava inmers en el càlcul el qui em cobrava diu:

$1.04+1.04+1.24 = 3.32$. Aquesta era fàcil, ehh!! Que necessitaves una calculadora tu??

I el meu company hi afegeix que el que necessitava és el meu ordinador, que sense ordinador no sé calcular xD

I automàticament vaig començar a pensar com calcular-ho amb l'ajut d'un ordinador però sense haver de fer cap suma, clar...

La resposta?? Monte Carlo. D'aquest tema ja en van escriure quatre ratlles introductòries a Cocociència. El que presentaré aquí és com vaig pensar que es podia fer la suma.

La idea dels mètodes de Monte Carlo és fer ús de nombres aleatòris per a calcular integrals. En cas de la suma es tractava de fer el càlcul d'una superfície. A l'eix x assignem un intèrval de longitud 1 a cada un dels producte mentre que a l'eix y hi assignem el preu per unitat tal i com es pot veure a la figura inferior.

És fàcil veure que si sumem les àrees vermella i verda obtindrem el valor total de la compra. Però això no ens resol el problema ja que per a fer-ho hauríem de sumar. La resposta es troba en el mètode de Monte Carlo.

La idea és molt simple. Agafem una parella de nombres aleatòris, $x$ i $y$. El primer que pugui prendre valors compresos entre el 0 i el 3, $x\in(0,3)$ i el segon que els prengui entre el 0 i el 2, $y\in(0,2)$. Amb aquesta parella de números podem pintar el punt amb coordenades $(x, y)$.

El que fem a continuació és veure si el punt $(x, y)$ ha caigut dins de la zona pintada o si ha caigut a fora. Si es troba dins de la zona pintada ens ho apuntem i comptem que hem encertat una vegada més.

Després de fer molts llençaments de números aleatòris haurem comptat quantes vegades hem encertat a la zona pintada, $n_{hit}$ i quants llençaments hem fet en total, $n_{tot}$.

Amb això podem saber quina fracció de llençaments, $P$, ha estat certera $P=n_{hit}/n_{tot}$.  Ara haurem de multiplicar aquesta fracció per l'àrea total que estàvem explorant, es a dir, la superfície, $S$, que ens era accessible en cada un dels llençaments. $S=2 \cdot 3$

Ara doncs serà fàcil de calcular el valor de la suma. Aquest serà $6 P$.

El resultat, obtingut en 20 càlculs amb 50.000 tirades aleatòries és el següent:

 

I si el que representem és la mitjana sobre 20 simulacions independents obtenim el següent:

Aquestes corbes són "uniques" de cada un dels càlculs però a l'infinit haurien de tendir al valor de la integral que és, en aquest cas, 3.32

Com es pot veuure l'error estadístic després de tota aquesta història és d'un cèntim d'€, de l'1%. No està malament!!

Sèries: El problema de Basilea

Ara feia molt temps que no entrava al blog i escrivia alguna cosa. No és per falta de temes ja que sempre en tinc al cap, l'únic que passa és que em fa una mandra espantosa tornar a escriure. De mica en mica, però, aniré tornant. Ara per tornar he vist que el LaTeX tornava a fallar i l'he arreglat i en aquest article ficaré 4 números i sumes infinites per posar-lo a prova. Benvinguts del nou al meu bloc els qui seguiu per aquest món.

Quan un aplica les matemàtiques a la vida real en moltes ocasions es troba amb les sèries. Una sèrie és una suma de termes d'una successió. Algunes d'elles tenen la seva història. Fa temps vaig parlar de la sèrie aritmètica que la va sumar Gauss i avui toca el problema de Basilea que va ser resolt per Euler. Euler, com molts altres matemàtics després seu, va ser dels que podrien haver fet seva la frase del gran Buzz-Lightyear.

El problema de Basilea és fer la suma següent:

\[\sum_{n=1}^{\infty}\frac{1}{n^2}=1+\frac{1}{4}+\frac{1}{9}+\frac{1}{16}+...\]

A diferència de la suma aritmètica aquesta és completament impossible de fer a mà però sí que ens hi podem anar acostant de mica en mica. La idea per a fer la suma és buscar-la en un lloc conegut.

Per a fer-ho recorrem a les sèries de Taylor (que ens ajuden en infinitat de situacions). El desenvolupament en sèrie de Taylor del sinus és:

\[sin(x)=x-\frac{x^3}{6}+\frac{x^5}{120}-...=x(1-\frac{x^2}{6}+\frac{x^4}{120}-...)\]

I si dividim per x obtenim:

\[\frac{sin(x)}{x}=1-\frac{x^2}{6}+\frac{x^4}{120}-...\]

D'altra banda podem construir la funció anterior emprant un mètode més còmode. Els punts en què de la funció $sin(x)/x$ s'anul·la en els mateixos punts que la funció $\sin(x)$ i són $x=\pm\pi$, $x=\pm2\pi$... De fet té infinits punts en què s'anul·la, tots els què satisfan $x=n\pi$. Per tant, podem construir un polinomi que representi la nostra funció de la manera següent:

\[\frac{sin(x)}{x}=(1-\frac{x}{\pi})(1+\frac{x}{\pi})(1-\frac{x}{2\pi})(1+\frac{x}{\pi})(1-\frac{x}{3\pi})(1+\frac{x}{3\pi})...\]

Veiem que la funció s'anul·la en tots els punts que havíem mencionat abans $x=n\pi$. Ara podem aplicar la relació $(a+b)(a-b)=a^2-b^2$ i obtenim una expressió amb què ens serà més còmode de treballar.

\[\frac{sin(x)}{x}=(1-\frac{x^2}{\pi^2})(1-\frac{x^2}{2^2\pi^2})(1-\frac{x^2}{3^2\pi^2})...\]

Ara només hem de començar a multiplicar com a autòmats (hi ha altres expressions més freqüents però de moment els autòmats no se senten ofesos):

\[\frac{sin(x)}{x}=1-(\frac{x^2}{\pi^2}+\frac{x^2}{2^2\pi^2}+\frac{x^2}{3^2\pi^2}+...)+(\frac{x^2}{\pi^2}\frac{x^2}{2^2\pi^2}+...)-...=\]

\[=1-\frac{x^2}{\pi^2}(1+\frac{1}{2^2}+\frac{1}{3^2}+...)+\frac{x^4}{\pi^4}(\frac{1}{1^2}\frac{1}{2^2}+\frac{1}{1^2}\frac{1^2}{3^2}+...)-...\]

On el primer terme és el producte de tots els uns. El segon terme és multiplicar cada terme amb $x^2$ amb els uns, podem veure que anirem tenint termes de la forma $\frac{x^2}{n^2\pi^2}$. El tercer terme és el resultat de multiplicar dos termes amb $x^2$ diferents pels uns.. Fent-ho així acabaríem fent el producte de tots amb tots però ja els tindríem ordenats en potències parelles de $x$. De moment ens centrem en el segon terme: el que va amb $x^2$, és fàcil de veure que és la sèrie que preteníem sumar!!

\[\frac{1}{1^2}+\frac{1}{2^2}+\frac{1}{3^2}+\frac{1}{4^2}=\sum_{n=1}^{\infty}\frac{1}{n^2}\]

Ara doncs ens queda comparar els dos coeficients que acompanyen el terme $x^2$ dels dos desenvolupaments en sèrie.

\[-\frac{1}{6}=-\frac{1}{\pi^2}\sum_{n=1}^{\infty}\frac{1}{n^2}\longrightarrow\sum_{n=1}^{\infty}\frac{1}{n^2}=\frac{\pi^2}{6}\]

Amb aquest mètode, una mica indirecte, hem aconseguit fer una suma infinita sense haver hagut de sumar cap número. Amb accions com aquesta els grans matemàtics passen a la història.

Ara podem aprofitar la idea per a fer el càlcul amb $1/n^4$. Per a fer-ho jo he hagut de fer una resta cosa que havia de passar, no tinc el nivell de l'Euler.

El terme amb $x^4$ el podem expandir una mica més i obtenim:

\[\frac{x^2}{\pi^2}\frac{x^2}{2^2\pi^2}+\frac{x^2}{\pi^2}\frac{x^2}{3^2\pi^2}+...+\frac{x^2}{2^2\pi^2}\frac{x^2}{3^2\pi^2}+\frac{x^2}{2\pi^2}\frac{x^2}{4^2\pi^2}+...\]

\[...+\frac{x^2}{3\pi^2}\frac{x^2}{4^2\pi^2}+\frac{x^2}{3\pi^2}\frac{x^2}{5^2\pi^2}...=\]

\[=\frac{x^4}{\pi^4}(\frac{1}{1^2}(\frac{1}{2^2}+\frac{1}{3^2}+...)+\frac{1}{2^2}\left(\frac{1}{3^2}+\frac{1}{4^2}+...)+\frac{1}{3^2}(\frac{1}{4^2}+...)\right)\]

Podem veure com s'ordenen. Veiem que cada terme de la suma consta de dos termes, el primer recorre tots els valors possibles de n mentre que el segon recorre només els k termes que són més grans que n. En format compacte això ho podríem escriure de la manera següent:

\[\sum_{n=1}^{\infty}\frac{1}{n^2}\sum_{k>n}^{\infty}\frac{1}{k^2}=\sum_{n=1}^{\infty}\sum_{k>n}^{\infty}\frac{1}{n^2k^2}=\frac{\pi^2}{120}\]

On l'últim resultat és el de comparar amb el terme d'ordre $x^4$ del desenvolupament en sèrie. Ara és quan hem de mirar de sumar això d'alguna manera i mirar de fer servir el resultat anterior (El del problema de Basilea). Si miréssim de sumar per a tots els $n$ i $k$ sense restriccions estaríem ficant termes de més, concretament estaríem afegint tots els termes amb k més petit que n (terme igual que el $k>n$ (només hem de permutar els noms a n i k i es veu de seguida)) i un terme amb $n=k$.

\[\sum_{n=1}^{\infty}\sum_{k=1}^{\infty}\frac{1}{n^2k^2}=2\sum_{n=1}^{\infty}\sum_{k>n}\frac{1}{n^2k^2}+\sum_{n=1}^{\infty}\frac{1}{n^4}\]

I ara sí que és fàcil de fer la primera suma, la que no té limits sobre n ni k.

\[\sum_{n=1}^{\infty}\frac{1}{n^2}\sum_{k=1}^{\infty}\frac{1}{k^2}=\left(\sum_{n=1}^{\infty}\frac{1}{n^2}\right)^2=\left(\frac{\pi^2}{6}\right)^2=\frac{\pi^4}{36}\]

Ara ja només queda substituir valors i anirem a petar a una equació de primer grau:

\[\frac{\pi^4}{36}=2\frac{\pi^4}{120}+\sum_{n=1}^{\infty}\frac{1}{n^4}\]

\[\sum_{n=1}^{\infty}\frac{1}{n^4}=\frac{\pi^4}{36}-\frac{2\pi^4}{120}=\frac{10\pi^4}{360}-\frac{6\pi^4}{360}=\frac{\pi^4}{90}\]

Podem comprovar que és el resultat que s'espera... Fer els següents és senzill havent fet ja el cas $1/n^4$ però és més pesat, s'ha de reconèixer.

L'IFS i la falguera de Barnsley

Fa un temps en aquest blog vaig escriure sobre el joc del caos i petites variacions. Aquests jocs es basen en el que a dia d'avui es coneix com el Teorema del Collage. Si llegim el que en diu la Wikipedia no s'entén pràcticament res (és el que passa amb molts teoremes).

Si mirem d'entendre què passava amb el joc del caos veiem que tenim 3 punts especials (cadascun dels vèrtex) i que estem fent una reducció de longitud des d'un punt qualsevol a un dels vèrtex. D'haver-hi només un dels vèrtex el punt cauria, inevitablement, cap al vèrtex: el vèrtex és un atractor. En ficar els 3 vèrtex junts el punt no acaba de caure mai en cap d'ells (tot i que s'hi acosta i es manté acotat) el que diu el teorema del collage és que l'atractor total és la unió dels atractors independents. Costa d'entendre però mirant una estoneta el triangle de Sierpinski un se'n fa una idea més o menys encertada.

El matemàtic Michael Barnsley al llibre fractals everywhere ens proposa un exemple combinant 4 transformacions afins (una mica més complicades que les del joc del caos). Esquemàticament tenen aquesta pinta:

La primera, representada en verd és una reducció d'escala de la coordenada y, amb una probabilitat de l'1% de ser aplicada:

La segona, representada per un rectangle blau s'aplica un 85% de les vegades:

La tercera, representada per un rectangle vermell s'aplica en un 7% de les vegades:

La quarta, representada per un rectangle blau fosc s'aplica en un 7%:

Si fem que l'ordinador ho dibuixi obtenim, efectivament, la falguera de Barnsley. A continuació veurem diferents figures amb 10, 50, 100, 500, 1000, 10000 i 50000 iteracions.

Fes el teu fractal de Newton

Les imatges dels fractals de Newton de l'entrada passada els vaig fer amb el Mathematica. Ara bé, hi ha mètodes alternatius que són més ràpids, i còmodes de fer anar. L'únic inconvenient que té és que s'ha d'introduir la derivada a mà.

El mètode alternatiu implica el programa Ultra Fractal 5. Aquest és un programa que cal comprar o intal·lar-ne una versió de prova. La meva recomanació és que instal·leu la versió de prova ja que quan acaba el termini es pot seguir utilitzant. L'únic problema és que té inhabilitades algunes opcions de renderització i que deixa marques d'aigua a les exportacions d'imatges. Tanmateix és el programa ideal per jugar amb el pla complex.

Si us decidiu per intal·lar el programa podreu navegar pels conjunts de Julia i Mandelbrot i crear-vos els conjunts de Newton que vulgueu. I no només navegar ja que podreu jugar amb les opcions de coloració dels fractals, capes de colors, etc. Amb ell surt la vena artística de la gent de ciències!!

Acabada la publicitat (ara és quan aquella gent m'hauria de donar una llicència gratuïta pel programa) toca començar a fer fractals.

Si teniu ja el programa (i després d'obrir-lo i explorar el conjunt de Mandelbrot) es pot provar de fer fractals de Newton. Per a fer-ho només cal seguir les instruccions següents.

1.- Aneu a File, New, Fractal Formula File.

2.- Entreu a Insert, New Formula..

3.- Introduïu el nom de l'arxiu "Newton", per exemple. Després d'accepta apareix un text base amb comentaris. Ara només caldrà donar algunes instruccions bàsiques.

4.- A l'apartat init: s'ha d'introduir la intrucció z=#pixel. Això dóna a la variable z el valor del píxel en què es troba (punt inicial pel mètode de Newton).

5.- A l'apartat loop: s'ha d'introduir l'expressió del mètode de Newton per exemple z=z-(z^3+2*z^2+z+1)/(3*z^2+4*z+1). Això generarà el fractal de Newton corresponent a la funció $f(z)=z^3+2*z^2+z+1$.

6.- A l'apartat bailout: s'hi ha d'introduir la condició que atura el loop:. En el nostre cas volem que deixi de calcular quan f(z) sigui proper a zero. |x^3+2*z^2+z+1|>0.005 farà que el mètode deixi d'aplicar-se quan f(z) sigui més petita que 0.005.

7.- Ara només cal guardar la fórmula File, Save as, Formulas, My Formulas.

8.- Per obrir el fractal File, New, Fractal, My Formulas, Newton.ufm.

9.- I a explorar!

En un dels comentaris de l'article anterior en Sheldon va preguntar cap a on convergia cada punt ja que els dibuixos els vaig fer en funció de la velocitat de convergència. He estat explorant l'Ultra Fractal i aquí ve el resultat. Les zones estan pintades en funció de la solució a la que convergeixen.

Per la funció $f(z)=z^3+2z^2+z+1$ teníem clarament 3 zones de convergència i una frontera fractal, per tant obtenim aquestes figures:

Per $f(z)=z^3+2z^2+\sqrt{z}+1$:

I per $f(z)=z^3+2z^2+ln(z)+1$:

I ara un dibuix fet amb $f(z)=z^3+2z^2+tan(z)+1$:

Fractals i mètode de Newton

En l'entrada anterior vaig parlar del mètode de Newton per tal de trobar zeros de funcions i resoldre equacions. Com a exemple vaig trobar la solució real de $x^3+2x^2+x+1=0$. A més a més, vaig comentar que n'hi havia dues més i que eren complexes. Per tant, el més raonable és ficar-se a treballar amb els nombres complexos.

Per tenir una idea aproximada d'on es troben les solucions he anat avaluant la funció polinòmica en diferents punts del pla complex (propers a l'origen) i he pintat cada punt depenent del valor del seu mòdul. Els punts més foscos es corresponen als valors més alts.

Els tres zeros es corresponen a les zones blanques i són $x_1=-1.755$ (a l'esquerra), $x_2=-0.123+0.745i$ (a la part superior) i $x_3=-0.123-0.745i$ (a la part inferior). Fins aquí l'únic que s'ha hagut de fer és fer ús del mètode de Newton. El mètode de Newton necessita un punt d'inici, cap a quina solució anem depenent d'aquest punt? O són tots els punts igual de bons??

El que he fet a continuació és anar aplicant el mètode de Newton a diferents punts del pla complex i pintar-los en funció del nombre d'iteracions que necessiten per convergir a una solució.

Pintades de colors més foscos (lila) es mostren les zones per les quals la convergència és més ràpida. Com era d'esperar el pla queda clarament dividit en 3 regions els punts de les quals convergeixen cap a cada un dels zeros. El més interessant, però, és la frontera! Si ampliem una mica la imatge a prop de la frontera que separa les arrels imaginàries ens trobem amb formes molt curioses:

Podem apreciar l'autosimilitud que caracteritza moltes estructures fractals. Aquestes figures es coneixen com els fractals de Newton. Ara doncs, l'únic que fa falta és trobar una equació i provar de resoldre-la.

L'equació $sin(x)+x=0$ només té un zero real ($x=0$). Ara bé el sinus d'un nombre complex pot prendre valors més grans que 1 (cosa que no passa amb els reals). Si es representa el valor de la funció per a diferents punts del pla podem veure el zero real i els zeros complexos (sempre aparellats amb el seu complex conjugat):

Els primers 4 zeros complexos es troben a $x_1=4.212+2.251i$, $x_2= 4.212-2.251i$, $x_3=-4.212+2.251i$ i $x_4=-4.212-2.251i$. Un cop localitzats els zeros pintem els punts en funció de la convergència a les arrels i ampliem algunes zones.

Si canviem la funció per $f(x)=sin(x)+x^2$ obtenim altres estructures:

Si ens fixem amb funcions més complicades com les arrels quadrades i els logaritmes podem trobar-nos amb estructures més problemàtiques. Per exemple, si a l'equació inicial de tercer grau substituïm la $x$ per $\sqrt{x}$ obtenim una nova funció $f(x)=x^3+2x^2+\sqrt{x}+1$ que porta a un fractal de Newton molt diferent.

I fent $f(x)=x^3+2x^2+ln(x)+1$ obtenim aquestes figures:

Ara només és qüestió de trobar funcions interessants i anar-les representant. I encara hi ha qui diu que les matemàtiques són avorrides!! És una llàstima que les matemàtiques les pugui gaudir tan poca gent...

El mètode de Newton

Un dels capítols importants de les matemàtiques als instituts és la resolució d'equacions. Allà s'ensenyen fórmules per poder arribar a resoldre equacions de segon grau del tipus $ax^2+bx+c=0$.

Si un es fica a inventar-se equacions veu que de seguida passen a ser més complicades de fer. Per exemple, per una equació de tercer grau la fórmula ja és espantosa. Com que recordar aquests desenvolupaments (o saber-los deduir) és feixuc els matemàtics han buscat i trobat camins alternatius.

El que ens interessa és trobar el que s'anomenen zeros d'una funció. Es a dir els valors de $x$ pels que la funció s'anul·la ($f(x)=0$) aquesta és la nostra equació). Per fer això hi ha mètodes força elementals (anar provant valors (més o menys com la loteria de Nadal)) i d'altres de més refinats com el de Newton-Raphson (també anomenat de la tangent).

El mètode de Newton necessita una funció $f(x)$ (que descriu una corba en el pla) i un valor qualsevol $x_1$. El que hem de fer és trobar la recta tangent a la corba en $x_1$ i mirar on talla l'eix de les x: aquest punt de tall serà el valor $x_2$. Amb el segon valor de x hem de fer el que hem fet amb $x_1$ per tal de trobar punt $x_3$. I així successivament.

L'expressió que es fa servir per calcular el punt de tall amb l'eix de les x és la següent:

$x_{n+1}=x_{n}-\frac{f(x_n)}{f\prime(x_n)}$.

On $f'(x_n)$ és la derivada de la funció en el punt $x_n$ (la derivada d'una funció en un punt dóna el pendent de la recta tangent a la funció en aquell punt, es a dir, com de ràpid creix).

En l'animació es mostra com donada una corba se'n troba el zero. La gràcia del mètode és que és ràpid i el podem fer servir en moltes situacions.

Si volem resoldre una equació de tercer grau com la següent $x^3+2x^2+x+1=0$ el primer que hem de fer és construir la funció $f(x)$ i la seva derivada:

$f(x)=x^3+2x^2+x+1$

$f\prime(x)=3x^2+4x+1$

Si representem la funció podem veure que el zero es troba al voltant de -1.8.

Començarem amb $x_1=2$ per veure com funciona el mètode i anirem avaluant la funció després de cada pas.

$x_1=2$

$f(x_1)=19$

$x_2=1.09523$

$f(x_2)= 5.80805$

$x_3=0.44842$

$f(x_3)=1.94075$

$x_4=-0.12290$

$f(x_4)= 0.90545$

$x_5=-1.75814$

$f(x_5)= -0.0105$

$x_6=-1.75489$

$f(x_6)= 0.00004$

Com podem veure amb pocs passos hem obtingut un nombre que fa f(x) molt propera a 0. El punt $x_6$ és una bona aproximació de la solució de l'equació.

Tanmateix, les equacions de tercer grau tenen sempre 3 solucions tot i que poden no ser entre els nombres reals. La que tenim nosaltres només presenta un punt de tall amb l'eix x i, per tant, només té una solució real (les altres dues es troben en els nombres complexos). Per trobar les dues restants també es pot fer ús del mètode de Newton. Depenent de quin sigui el punt inicial anirem cap a zeros diferents!

Les preguntes que un sempre s'ha de fer amb els mètodes numèrics són les següents:

1.- Realment dóna el valor que volem??

2.- En cas de donar el valor correcte quan triga a arribar-hi?

Per tal de respondre aquestes preguntes ara que disposem d'ordinadors la resposta és òbvia. Agafem una regió i anem provant diferents nombres i ens anotem els passos que fan falta per obtenir un valor determinat. Però això ja és del proper article.