Monthly Archives: June 2010

by

Ante la duda, aplicar Fourier

Introducción

Un clásico problema para explicar el poder del principio de superposición y de la simetría consiste en calcular la resistencia entre los puntos A y B de la siguiente red infinita de resistores:

Red infinita de resistores unitarios con dos nodos marcados.

La forma más simple de resolverlo es asumir que inyectamos una corriente de 1 por A, dejando desconectados los otros nodos. Por simetría, las corrientes que partan del nodo serán todas iguales a 1/4:

Se inyecta una corriente unitaria en el nodo A de la red anterior.

Por linealidad y simetría, si repetimos el procedimiento inyectando una corriente de -1 por el nodo B, las corrientes que partan del mismo serán iguales a -1/4. Aplicando el principio de superposición (posible por la linealidad de la Ley de Ohm), llegamos a las siguientes corrientes:

Superposición de los efectos de inyectar una corriente unitaria por el nodo A y extraer una corriente unitaria por el nodo B.

Como esto equivale a inyectar una corriente de 1 por A y extraerla por B, la diferencia de potencial entre ambos nodos será idéntica a la resistencia entre estos. La corriente que circula por el resistor unitario que los une es 1/2 y, por lo tanto, la resistencia entre A y B es 1/2.

Una variante de este problema (considerablemente más compleja) que ha aparecido en múltiples lugares (incluso XKCD) es la siguiente:

Red de resistores unitarios con los nodos A y B separados por un 'salto de caballo'.

En lo que resta del post veremos como obtener una solución numérica y, finalmente, una analítica a este problema.

Una solución numérica

Un nodo cualquiera de la red de resistores debe satisfacer la ley de Kirchhoff de conservación de corrientes. Como las corrientes a través de cada resistor unitario son iguales a la diferencia de tensión, tenemos:

(V_{i\,j} - V_{i-1\,j}) + (V_{i\,j} - V_{i+1\,j}) + (V_{i\,j} - V_{i\,j-1}) + (V_{i\,j} - V_{i\,j+1}) = 0

4V_{i\,j} - V_{i-1\,j} - V_{i+1\,j} - V_{i\,j-1} - V_{i\,j+1} = 0.

Estas ecuaciones se aplicarán a todos los nodos a excepción de los que estén conectados a fuentes externas. Si reemplazamos la red infinita de resistores por una red finita con N2 nodos, la aplicación de la ecuación de nodos nos dejará N2 – 2 ecuaciones, ya que la ecuación de nodos no se aplicará donde conectemos la fuente. Esta cantidad de ecuaciones es suficiente para determinar todas las tensiones, ya que las tensiones de dos nodos son forzadas por la fuente.

Utilizando una topología toroidal por simplicidad y aplicando Gauss Seidel + SOR con ω = 1.5 llegamos al siguiente código (ver resultados en codepad):

TOL = 1e-6
K = 1.5

def solve(n):
    assert n >= 4
    v = [[0.5 for y in range(n)] for x in range(n)]
    delta = 2 * TOL
    FIXED_VALUES = {
        (0, 0): 1,
        (2, 1): 0
    }
    while delta > TOL:
        delta = 0
        for x in range(n):
            for y in range(n):
                if (x, y) in FIXED_VALUES:
                    v[x][y] = FIXED_VALUES[(x, y)]
                    continue
                vo = v[x][y]
                vn = (v[x - 1][y] + v[(x + 1) % n][y] +
                      v[x][y - 1] + v[x][(y + 1) % n]) / 4.0
                v[x][y] = vn * K + vo * (1 - K)
                delta = max(abs(vn - vo), delta)
    return 1.0 / (4 * v[0][0] - v[-1][0] - v[1][0] - v[0][-1] - v[0][1])

if __name__ == '__main__':
    from time import time
    for n in (8, 16, 32, 64, 128):
        start = time()
        print 'N = %d - R = %f - Time: %f s' % (n, solve(n), time() - start)

Ejecutando el código completo (la versión de codepad excluye N = 64 y N = 128 por las limitaciones al tiempo de ejecución) obtenemos la siguiente salida:

N = 8 - R = 0.734997 - Time: 0.000000 s
N = 16 - R = 0.763460 - Time: 0.094000 s
N = 32 - R = 0.770537 - Time: 0.781000 s
N = 64 - R = 0.772408 - Time: 2.532000 s
N = 128 - R = 0.773033 - Time: 23.687000 s

Se observa una convergencia relativamente rápida a un valor cercano a 0.773 y un comportamiento bastante extraño del tiempo de ejecución.

Una solución analítica

Si llamamos Vi j a la tensión del nodo m, n e Im n a la corriente inyectada en el mismo, la conservación de la corriente nos lleva a lo siguiente:

I_{m\,n} = 4V_{m\,n} - V_{m-1\,n} - V_{m+1\,n} - V_{m\,n-1} - V_{m\,n+1}.

Las corrientes inyectadas pueden representarse como

I_{m\,n} = \delta_{m\,n} - \delta_{m-M\,n-N},

donde las “deltas” inyecta una corriente unitaria en el nodo 0, 0 y la extraen por el nodo M, N. Esta ecuación en diferencias no tiene una solución única por esencialmente el mismo motivo que el potencial de una carga puntual no tiene una única solución: el núcleo del laplaciano es no trivial. Pero hay una solución única si exigimos que el potencial se vaya a cero cuando m y n se van a infinito.

Un método muy utilizado para resolver ecuaciones diferenciales es aplicar una transformación inversible que convierta a la ecuación diferencial en una simple ecuación algebraica. Si esta última puede resolverse, aplicando la transformación inversa obtendremos una solución en el dominio original.

Si aplicamos la “versión adecuada” de la transformada de Fourier en dos dimensiones,

\displaystyle V(p,q) = \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty} V_{m\,n} e^{-i m p} e^{-i n q},

podemos ver que la transformada de Vm-1 n es

\displaystyle \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty} V_{m-1\,n} e^{i m p} e^{i n q} = \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty} V_{m\,n} e^{-i (m + 1) p} e^{-i n q}

\displaystyle = e^{-i p} \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty} V_{m\,n} e^{-i m p} e^{-i n q}

\displaystyle = e^{-i p} V(p, q).

Utilizando esto, llegamos a la siguiente ecuación algebraica sobre V(p, q):

\displaystyle 1 - e^{-i M p}e^{-i N q} = V(p, q) (4 - e^{-i p} - e^{i p} - e^{-i q} - e^{i q}).

Resolviéndola:

\displaystyle V(p, q) = \frac{1 - e^{-i M p}e^{-i N q}}{4 - 2 \cos p - 2 \cos q}.

Antitransformando:

\displaystyle V_{m\,n} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - e^{-i M p}e^{-i N q}}{4 - 2 \cos p - 2 \cos q}e^{i m p}e^{i n q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - e^{-i M p}e^{-i N q}}{4 - 2 \cos p - 2 \cos q}e^{i m p}e^{i n q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{e^{i m p}e^{i n q} - e^{i (m - M) p}e^{i (n - N) q}}{4 - 2 \cos p - 2 \cos q}.

Pero a nosotros nos interesa la resistencia entre los nodos 0, 0 y M, N, igual a la diferencia de tensión entre ambos nodos (ya que circula entre ellos una corriente unitaria). Por lo tanto,

\displaystyle R_{M\,N} = V_{0\,0} - V_{M\,N}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - e^{-i M p}e^{-i N q} - e^{i M p}e^{i N q} + 1}{4 - 2 \cos p - 2 \cos q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - e^{-i (M p + N q)} - e^{i (M p + N q)} + 1}{4 - 2 \cos p - 2 \cos q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{2 - 2\cos(M p + N q)}{4 - 2 \cos p - 2 \cos q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos(M p + N q)}{2 - \cos p - \cos q}

Ahora “solo” resta resolver esta integral para cuando M = 2 y N = 1:

\displaystyle R_{2\,1} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos(2 p + q)}{2 - \cos p - \cos q}.

Evaluándola numéricamente obtenemos 0.7732, un valor que coincide en 3 cifras significativas con el obtenido previamente.

Solución analítica de la integral

ADVERTENCIA: gran cantidad de fórmulas que no aportan mucho…

Podemos empezar buscando una solución via Wolfram Alpha, pero no funciona. Entonces comenzamos expandiendo el numerador:

\displaystyle R_{2\,1} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos(2 p + q)}{2 - \cos p - \cos q}

\displaystyle = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos 2p \cos q + \sin 2p \sin q}{2 - \cos p - \cos q}

Trabajando con el primer término,

\displaystyle A = \frac{1}{2 - \cos p}

\displaystyle \int_0^{2\pi} dq \frac{1}{2 - \cos p - \cos q} = \frac{1}{A^{-1}}\int_0^{2\pi} dq \frac{1}{1 - A\cos q},

llegamos a una forma en tabla (Tabla 64, Formula 12):

\displaystyle = 2A\int_0^{\pi} dq \frac{\cos 0\cdot q}{1 + (-A)\cos q}

Aplicando la fórmula y simplificando:

\displaystyle = 2A\frac{\pi}{\sqrt{1-A^2}}\left(\frac{\sqrt{1-A^2}-1}{-A}\right)^0

\displaystyle = \frac{2}{2 - \cos p}\frac{\pi}{\sqrt{\frac{(2 - \cos p)^2-1}{(2 - \cos p)^2}}}

\displaystyle = \frac{2}{2 - \cos p}\frac{\pi(2 - \cos p)}{\sqrt{(2 - \cos p)^2-1}}

\displaystyle = \frac{2\pi}{\sqrt{(2 - \cos p)^2-1}}

Del mismo modo podemos proceder con el segundo término:

\displaystyle \int_0^{2\pi} dq \frac{\cos 2p \cos q}{2 - \cos p - \cos q} = \frac{\cos 2p}{A^{-1}}\int_0^{2\pi} dq \frac{\cos q}{1 - A\cos q}

\displaystyle = A \cos 2p\int_0^{2\pi} dq \frac{\cos q}{1 - A\cos q}

\displaystyle = 2A\cos 2p\int_0^{\pi} dq \frac{\cos 1\cdot q}{1 + (-A)\cos q}

\displaystyle = 2A\cos 2p\frac{\pi}{\sqrt{1-A^2}}\left(\frac{\sqrt{1-A^2}-1}{-A}\right)^1

\displaystyle = -2A\cos 2p\frac{\pi}{\sqrt{1-A^2}}\frac{\sqrt{1-A^2}-1}{A}

\displaystyle = -2\frac{1}{2 - \cos p}\cos 2p\frac{\pi}{\sqrt{\frac{(2 - \cos p)^2-1}{(2 - \cos p)^2}}}\frac{\sqrt{\frac{(2 - \cos p)^2-1}{(2 - \cos p)^2}}-1}{\frac{1}{2 - \cos p}}

\displaystyle = -\frac{2\pi\cos 2p}{(2 - \cos p)\frac{\sqrt{(2 - \cos p)^2-1}}{2 - \cos p}}\frac{\frac{\sqrt{(2 - \cos p)^2-1}}{2 - \cos p}-1}{\frac{1}{2 - \cos p}}

\displaystyle = -\frac{2\pi\cos 2p}{\sqrt{(2 - \cos p)^2-1}}\left(\sqrt{(2 - \cos p)^2-1} - (2 - \cos p)\right)

\displaystyle = - 2\pi\cos 2p + \frac{2\pi\cos 2p(2 - \cos p)}{\sqrt{(2 - \cos p)^2-1}}

La integral del último término se anula, como puede verse a continuación:

\displaystyle \int_0^{2\pi} dq \frac{\sin 2p \sin q}{2 - \cos p - \cos q} = \frac{1}{A^{-1}}\int_0^{2\pi} dq \frac{\sin 2p \sin q}{1 - A\cos q}

\displaystyle = A\sin 2p \int_0^{2\pi} dq \frac{\sin q}{1 - A\cos q}

\displaystyle = A\sin 2p \int_0^\pi dq \frac{\sin q}{1 - A\cos q} - A\sin 2p \int_{2\pi}^\pi dq \frac{\sin q}{1 - A\cos q}

\displaystyle = A\sin 2p \int_0^\pi dq \frac{\sin q}{1 - A\cos q} + A\sin 2p \int_0^\pi dr \frac{\sin (2\pi - r)}{1 - A\cos (2\pi - r)}

\displaystyle = A\sin 2p \int_0^\pi dq \frac{\sin q}{1 - A\cos q} + A\sin 2p \int_0^\pi dr \frac{-\sin r}{1 - A\cos r}

\displaystyle = 0

Mediante estas integraciones, la integral doble original

\displaystyle R_{2\,1} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \int_0^{2\pi} dq \frac{1 - \cos 2p \cos q + \sin 2p \sin q}{2 - \cos p - \cos q}

se reduce a una integral simple:

\displaystyle R_{2\,1} = \frac{1}{4\pi^2} \int_0^{2\pi} dp \frac{2\pi}{\sqrt{(2 - \cos p)^2-1}} - \frac{2\pi\cos 2p(2 - \cos p)}{\sqrt{(2 - \cos p)^2-1}} + 2\pi\cos 2p

\displaystyle = \frac{1}{2\pi} \int_0^{2\pi} dp \frac{1- 2\cos 2p + \cos p\cos 2p}{\sqrt{(2 - \cos p)^2-1}}.

Utilizando Wolfram Alpha:

\displaystyle R_{2\,1} = \frac{1}{2\pi} \int_0^{2\pi} dp \frac{1- 2\cos 2p + \cos p\cos 2p}{\sqrt{(2 - \cos p)^2-1}}

\displaystyle = \frac{1}{2\pi} (8-\pi)

\displaystyle = \frac{4}{\pi} - \frac{1}{2}

\displaystyle \approx 0.7732,

coincidiendo con el resultado aproximado anterior.

Conclusiones y reconocimientos

Esto vindica, como Fer bien conoce, la importancia de aplicar Fourier cuando uno se encuentra frente a un problema que no sabe como resolver. 😛

Soluciones esencialmente idénticas a esta fueron desarrolladas mucho antes por las fuentes citadas y por “frooha”, que desarrolló una solución muy similar a esta, entre otros. También existe al menos una solución por la via experimental, preferida por los electrónicos.

Advertisement