Lab 3 Solucin de sistema de ec lineales

Lab. 3: Solución de sistema de ec. lineales con Montecarlo / Metropolis Algoritmos Paralelos Glen Rodríguez

Aplicaciones n n Redes eléctricas Redes de tránsito Cadenas de Markov Solución de algunos tipos de ecuaciones diferenciales parciales

Problema n n Obtener vector x (tamaño nx 1) que cumple: Ax=b. Donde b es un vector conocido del mismo tamaño de x, y A es una matriz de tamaño nxn que cumple: ¡ ¡ A es NO SINGULAR Se puede escoger una matriz M no singular tal que MA=I-L, y que |λ(L)| < 1 para todos los eigenvalues λ de L (no singular también)

Algoritmo n n λ es un número real llamado eigenvalue si hay solución de Lu= λu El problema se puede ver como: ¡ n n x=Lx + Mb = Lx + f Si M=I, se puede demostrar que la sumatoria ΣLi converge a A-1 si i ∞ En los ‘ 50 s se demostró que un “random walk” puede aproximar esta sumatoria.

Algoritmo n x 0=0; xm=b+Lxm-1 donde xm=m-esima aproximación a x

l ij

Pseudocódigo 1. 2. 3. 4. 5. 6. 7. 8. 9. Calcular pij, vij. Se puede asumir que pi=constante (ej. : 0. 01). Entonces se puede hacer pij iguales en una fila i =(1 -pi)/numero de pij ≠ 0 en fila i=ti; vij=lij/ti. O se puede hacer pij= (1 -pi)* | lij | /Σ |lij| en fila i; vij= Σ |pij| en fila i / (1 -pi) For elem=1 to n x_sum=0; for k=1 to m v_prd=1; ; ult=entero entre 1 y n tal que lult, elem ≠ 0; x_est=bult/pult Escoger r= real aleatorio unif. entre 0 y 1 While r<=1 -pult Asociar r a un valor “escog” usando las pult, j v_prd= v_prd* vult, escog

Pseudocódigo 9. 10. 11. 12. 13. 14. 15. x_est=v_prd*bescog/pescog ult=escog end while x_sum= x_sum+x_est; end for x[elem]= x_sum/m; end for Que se puede paralelizar: cualquiera de los for. Mejoras: una vez conocido x[1], …x[q], usar esos valores para x[q+1], …, x[n]

Redes eléctricas n n Sólo necesitamos 3 leyes: Ohm: ΔV=IR entre los extremos de un conductor de resistencia R 1 ra de Kirchhoff: Σ Iin = Σ Iout en cada unión entre 2 o más conductores 2 da de Kirchhoff: Σ ΔV = ΣFuentes en un loop ó lazo cerrado.

Ejemplo

Redes de tráfico n Se usa solo la 1 ra Kirchhoff.

Resolviendo Ec. de Poisson con Montecarlo n Sea la ecuación de Poisson: n Con fuente sinusoidal : -2π2 sen(π x) sen(π y) Y cero en todas las fronteras En diferencias finitas: n n

Solución con “random walks” n n n De un punto (i, j) comenzar el “random walk”. La probabilidad de moverse a cualquiera de los 4 puntos vecinos es ¼. Generar un número aleatorio para escoger el vecino. Añadir g(x, y) en la nueva posición Repetir hasta llegar a una frontera Esto fue solo UNA “random walk” Después de N “random walks” el estimado de u(i, j) es:

Paralelización de esta solución n n Hacer el siguiente proceso para todos los puntos Comenzar en una frontera ¡ ¡ ¡ De afuera hacia adentro Fila por fila Actualizar (intercambiar) data en el límite entre procesadores Actualizar el walk dentro de los puntos de un procesador Si el walk sale de un procesador, informarle al otro procesador

Metrópolis n n 1. 2. 3. Se usa para integrales con funciones de peso. Genera un random walk pero guiado por una función de peso w(x). Ejemplo en 2 -D comenzando de un punto (xi, yi): Escoger δ, el step size Generar dos aleatorios R 1, R 2 uniformes en el rango [ -1, 1] El nuevo punto será 1. 2. 4. x. Ti+1 = xi + δR 1 y. Ti+1 = yi + δR 2 Evaluar el ratio de w en el punto actual vs. el punto anterior 1. r= w(x. Ti+1 , y. Ti+1) / w(xi, yi)

Metrópolis 5. Si r>1 aceptar el nuevo punto y regresar a (2) 1. 6. 7. xi+1 = x. Ti+1 ; yi+1 = y. Ti+1 Si r<1 aceptar el nuevo punto con probabilidad r. O sea, generar un aleatorio uniforme η en [0, 1] y aceptar el punto solo si r>η (hacer lo mismo que en el paso 5) Caso contrario, rechazar el nuevo punto y volver a (2) El step size no debe ser ni muy chico ni muy grande. Los puntos obtenidos se usan luego para integrar, como en Monte. Carlo (Constante*n/N). La paralelización es similar a la discutida para las integrales con Montecarlo