Una vez resuelta la ecuación de Laplace en coordenadas esféricas, por el método de separación de variables: $$\nabla ^2 V(r, \theta, \phi) = 0$$ se obtienen los siguientes modos normales: $$ V_{lm}(r, \theta, \phi) = \left( D_{lm} r^l + E_{lm} \frac{1}{r ^{-(l+1)}} \right) Y_l ^m(\theta, \phi) $$ Escritos también de la forma: $$ V_{lm}(r, \theta, \phi) = \left( D_{lm} r^l + E_{lm} \frac{1}{r ^{-(l+1)}} \right) e^{im\phi} P_l ^m (cos \theta). $$

Una vez evaluada las condiciones de frontera, es posible corroborar qué: $$ V(r, \theta, \phi) = \sum_ {l= 0} ^ \infty \sum_{-l} ^{l} \left( D_{lm} r^l \right) e^{im\phi} P_l ^m (cos \theta), ~~ l + m = 2k + 1 $$ Donde la solución divergente en el origen ha sido removida para garantizar que el potencial sea finito en el punto (0,0,0) contenido dentro del plano z = 0. La condición en la frontera de $z = 0 \rightarrow V(r, \theta, \phi) = 0$, restringe los posibles valores de $l$ y $m$ aquellos tales que $l + m = 2k + 1$.

Imponiendo la siguiente condición a la frontera, $$ \sum_ {l= 0} ^ \infty \sum_{-l} ^{l} \left( D_{lm} a^l \right) Y_l ^m(\theta, \phi) = f(\theta, \phi) $$ Donde: $$ f(\theta, \phi) = \left\{ \begin{array}{cc} V_0 sin{\frac{\phi}{2}} ~ cos^2 \theta ~ sin{\frac{\phi}{2}} & 0 \leq \theta \leq \frac{\pi}{2} \\ 0 & \frac{\pi}{2} \leq \theta \leq \pi \end{array} \right. $$ Multiplicando por $\hat{Y}_{l'} ^{m'}$: $$ \sum_ {l= 0} ^ \infty \sum_{-l} ^{l} \left( D_{lm} a^l \right) \hat{Y}_{l'} ^{m'} Y_l ^m = f(\theta, \phi) ~\hat{Y}_{l'} ^{m'}$$ Integrando: $$ \sum_ {l= 0} ^ \infty \sum_{-l} ^{l} \left( D_{lm} a^l \right) \int_0 ^\pi \int_0 ^{2\pi} \hat{Y}_{l'} ^{m'} Y_l ^m ~sin\theta d\theta d\phi= V_0 \int_0 ^{\frac{\pi}{2}} \int_0 ^{2\pi} sin{\frac{\phi}{2}} ~ cos^2 \theta ~\hat{Y}_{l'} ^{m'} ~ sin\theta d\theta d\phi + \int_\frac{\pi}{2} ^\pi \int_0 ^{2\pi} 0 ~sin\theta d\theta d\phi$$ Como: $$ \int_0 ^\pi \int_0 ^{2\pi}\hat{Y}_{l'} ^{m'} Y_l ^m ~sin\theta d\theta d\phi = \delta_{l, l'} \delta_{m, m'} $$

$$ D_{lm} a^l = V_0 \int_0 ^{\frac{\pi}{2}} \int_0 ^{2\pi} sin{\frac{\phi}{2}} ~ cos^2 \theta ~\hat{Y}_{l'} ^{m'} ~ sin\theta d\theta d\phi $$

Y siendo $\hat{Y}_l ^m (\theta, \phi) = N ~e^{-im\phi} P_l ^m (cos \theta)$, con $N = (-1)^m \sqrt{\frac{(2l + 1)(l-m)!}{4\pi (l+m)!}} $;

$$ D_{lm} a^l = V_0 N \int_0 ^{\frac{\pi}{2}} \int_0 ^{2\pi} sin{\frac{\phi}{2}} ~ cos^2 \theta ~e^{-im\phi} P_l ^m (cos \theta) ~ sin\theta d\theta d\phi $$

Separando las componentes en $\theta$ y $\phi$: $$ D_{lm} = \frac{V_0 N}{a^l} \left[ \int_0 ^{\frac{\pi}{2}} cos^2 \theta ~sin\theta ~P_l ^m (cos \theta) d\theta \right] \left[\int_0 ^{2\pi}sin{\frac{\phi}{2}} ~ e^{-im\phi} d\phi \right] $$

Con el cambio de variable: $x = cos\theta$, $- sin \theta d\theta = dx$. \begin{equation} D_{lm} = \frac{V_0 N}{a^l} \left[ -\int_1 ^0 x^2 ~P_l ^m (x) dx \right] \left[\int_0 ^{2\pi}sin{\frac{\phi}{2}} ~ e^{-im\phi} d\phi \right] \end{equation} Concentrándonos en la segunda integral, recordando que $sin \alpha = \frac{e^{i \alpha} - e ^{-i \alpha}}{2i} $; $$ \int sin{\frac{\phi}{2}} ~ e^{-im\phi} d\phi = \int \left( \frac{e^{i \frac{\phi}{2}} - e^{-i \frac{\phi}{2}} }{2i} \right) e ^{-i m \phi} d\phi $$

$$ = \frac{1}{2i} \int \left( e^{i\frac{\phi}{2}} e ^{-im\phi} - e^{-i \frac{\phi}{2}} e ^{-im\phi} \right)d\phi $$$$ = \frac{1}{2i} \int \left( e ^{i \phi \left( \frac{1}{2} - m \right) } - e ^{- i \phi \left( \frac{1}{2} + m \right)} \right) d\phi $$

Separando la integral y efectuando los cambios de variable pertinentes:

$$ = \left. \frac{1}{2} \left\{ \frac{1}{ \left( \frac{1}{2} + m \right) } e ^{- i \phi \left( \frac{1}{2} + m \right)} - \frac{1}{ \left( \frac{1}{2} - m \right) } e ^{i \phi \left( \frac{1}{2} - m \right)} \right\} \right|_0 ^{2\pi}$$$$ = \frac{1}{2} \left\{ \frac{1}{ \left( \frac{1}{2} + m \right) } e ^{- i \pi \left( 1 + 2m \right)} - \frac{1}{ \left( \frac{1}{2} - m \right) } e ^{i \pi \left( 1 - 2m \right)} \right\} - \frac{1}{2} \left\{ \frac{1}{ \left( \frac{1}{2} + m \right) } - \frac{1}{ \left( \frac{1}{2} - m \right) } \right\}$$

Usando la fórmula de Euler:

$$e ^{i \pi \left( 1 +2 m \right) } = cos ((2m +1) \pi) + i ~sin ((2m+1) \pi)$$

Sabemos que $m$ es un entero, por lo tanto $(2m+1)\pi$ es un múltiplo entero \emph{impar} de $\pi$ y $sin ((2m+1) \pi) = 0, \forall m$. Con el mismo argumento, $cos((2m+1)\pi) = -1$.

Análogamente: $$e ^{i \pi \left( 1 -2 m \right) } = cos ((1 - 2m) \pi) - i ~sin ((1 - 2m) \pi)$$ $$ = cos((2m-1) \pi) + i ~sin ((2m-1)\pi)) $$ Y por razones idénticas a las expuestas con anterioridad: $sin((2m-1)\pi) = 0$ y $cos((2m-1)\pi) = -1$, $\forall m$.

Hecho esto: $$ = \frac{1}{2} \left\{ \frac{1}{ \left( \frac{1}{2} + m \right) } (-1) - \frac{1}{ \left( \frac{1}{2} - m \right) } (-1) \right\} - \frac{1}{2} \left\{ \frac{1}{ \left( \frac{1}{2} + m \right) } - \frac{1}{ \left( \frac{1}{2} - m \right) } \right\}$$ $$ = - \frac{1}{2} \left\{ \frac{1}{ \left( \frac{1}{2} + m \right) } - \frac{1}{ \left( \frac{1}{2} - m \right) } \right\} - \frac{1}{2} \left\{ \frac{1}{ \left( \frac{1}{2} + m \right) } - \frac{1}{ \left( \frac{1}{2} - m \right) } \right\}$$ $$ = - \left\{ \frac{1}{ \left( \frac{1}{2} + m \right) } - \frac{1}{ \left( \frac{1}{2} - m \right) } \right\} $$

Y la expresión para los $D_{lm}$ se convierte en: $$ D_{lm} = \left\{ (-1)^m \sqrt{\frac{(2l + 1)(l-m)!}{4\pi (l+m)!}} \frac{V_0}{a^l} \left[\frac{1}{ \left( \frac{1}{2} - m \right) } - \frac{1}{ \left( \frac{1}{2} + m \right) } \right] \right\} \int_0 ^1 x^2 ~P_l ^m (x) dx $$

Definiendo $M_{lm}$: $$ M_{lm} = (-1)^m \sqrt{\frac{(2l + 1)(l-m)!}{4\pi (l+m)!}} \frac{V_0}{a^l} \left[\frac{1}{ \left( \frac{1}{2} - m \right) } - \frac{1}{ \left( \frac{1}{2} + m \right) } \right] $$

Es necesario notar que, para valores de $m = 0 $, los coeficientes $M_{lm} = 0$.


In [291]:
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm, colors
import scipy as spy
import scipy.integrate as integrate
from types import *
from mpl_toolkits.mplot3d import axes3d, Axes3D 
import matplotlib.cm as cm

%matplotlib inline

Introducimos los valores iniciales del problema. Por simplicidad, eligiremos una esfera unitaria y un potencial eléctrico de 5 volts.


In [292]:
a = 1.   #Radio de la Esfera
V0 = 5.  #Potencial eléctrico V0

Definimos una clase que permitirá efectuar operaciones con el grado y el orden de los polinomios con mayor eficiacia.


In [293]:
class Grado_Pol:
    
    def __init__ (self, l):
        if float(l).is_integer():                     # Comprueba que el grado sea un entero
            self.Grado = l
        else:
            print("El grado es invalido")
    
    def Orden (self):
        lista = []
        Emes = np.arange(-self.Grado, self.Grado +1)
        for m in Emes:
            if (self.Grado + m)%2 != 0:               # Elige unicamente los valores de m que satisfacen l+m = impar
                lista.append(m)
                
        return lista                                 # Devuelve una lista de valores posibles
    
    def No_Orden (self):
        return len(self.Orden())                      # Indica la cardinalidad de la lista anterior

A continuación, las funciones especiales que serán necesarias: Los polinomios asociados de legendre y los armonicos esferico, de los que únicamente será posible retomar la parte real.


In [294]:
def Legendre (x, grado, orden):
        return spy.special.lpmv(orden, grado, x)
def Ylm (r, th, phi, grado, orden):
    return spy.special.sph_harm(orden, grado, phi, th)

Definimos los coeficientes $D_{lm}$. La integral $\int_0 ^1 x^2 ~P_l ^m (x) dx $ será realizada con la libreria scypy.integrate.quad


In [295]:
def Clm (l, m):                                  # El coeficiente Clm es un auxiliar, permite conocer informacion util
    
    def Mlm (l, m):
        return ((-1.)**m)*np.sqrt(((2.*l + 1)*np.math.factorial(l-m))/(4*np.pi*np.math.factorial(l+m)))*(V0 / (a**l))*((1. / (0.5 - m)) - (1. / (0.5 + m)))
    integral = integrate.quad(Legendre, 1., 0., args=(l,m))
    
    return Mlm(l,m)*integral[0], integral[1], Mlm(l,m)
    # Clm[0] = el coeficiente buscado
    # Clm[1] = error del metodo quad de scipy
    # Clm[2] = Valor de Mlm, para ver si efectivamente es cero.

def Dlm (l,m):
    return Clm(l,m)[0]

Para la funcion potencial, es necesario realizar la "suma infinita" sobre los grados $l$ del polinomio. Por simplicidad, la variable de grado a sido renombrada como Rango, y puede manipularse para tomar una mejor aproximacion. Como primera propuesta, para ahorrar recursos de computo, ha sido designada con un valor predeterminado de 4.


In [296]:
def V(r, th, phi, Rango=4):
    V = 0.                           #Fuera del ciclo. Queremos sumar las todas las m's, sobre todas las i's.
    for i in range (1,Rango+1):
        l = Grado_Pol(i)
        if not l.Orden():
            print("El grado ", l.Grado, " no es valido")
        else:
            #print("Para el grado ", i, " hay ", l.No_Orden(), " ordenes validos.")
            for m in l.Orden():
                D = Dlm(i,m)
                if D < 1e-6 :
                    D = 0
                V = V + D*Ylm(r, th,phi, i, m)
                #print (D, "  l = ", i, " m = ", m, "  Error = ", Clm(i,m)[1], "Mlm = ", Clm(i,m)[2] )

                            
    return V

La funcion a graficar es, entonces, una funcion de 4 variables: tres coordenadas espaciales (en coordenadas esfericas), y una cuarta coordenada que representa el valor del potencial $V(r, \theta, \phi)$. Para visualizar una imagen, la variable $r$ será fijada como parámetro. La informacion del potencial sera incluida como un mapa de colores por medio de las librerias matplotlib.cm y matplotlib.colors, e incorporada a una plot_surface de la función $f(r, \theta, \phi)$ como una máscara de color, para cada punto.

Considerando únicamente la componente real de V y mapeando los puntos para Th y Ph obtenidos (un arreglo de N por N. De incluir la variable r sería de N x N x N):


In [297]:
"""
color_map= np.real(V(a,Th,Ph,5))
norm = colors.Normalize()
color = norm(color_map)

colorbar = cm.ScalarMappable(cmap = cm.jet)
colorbar.set_array(color_map)
"""


Out[297]:
'\ncolor_map= np.real(V(a,Th,Ph,5))\nnorm = colors.Normalize()\ncolor = norm(color_map)\n\ncolorbar = cm.ScalarMappable(cmap = cm.jet)\ncolorbar.set_array(color_map)\n'

Se hace necesario crear una malla de puntos para: $$ 0 \leq \theta \leq \pi/2 ~; ~ 0 \leq \phi \leq \pi$$


In [298]:
#Estos espacios son necesarios para la creacion de una malla de C x C puntos
C = 100

In [299]:
rr = np.linspace(0, a, C)
tt = np.linspace(0, (np.pi)/2. , C)
pp = np.linspace(0, 2*(np.pi), C)
RR,TT, PP = np.meshgrid(rr, tt, pp)        
Th,Ph = np.meshgrid( tt, pp)

Incluyendo el mapa de colores:


In [300]:
"""
** Esto tambien, necesito que corra para diferentes colores del radio

color_map= np.real(V(R,Th,Ph,5))
norm = colors.Normalize()
color = norm(color_map)

colorbar = cm.ScalarMappable(cmap = cm.jet)
colorbar.set_array(color_map)
"""


Out[300]:
'\n** Esto tambien, necesito que corra para diferentes colores del radio\n\ncolor_map= np.real(V(R,Th,Ph,5))\nnorm = colors.Normalize()\ncolor = norm(color_map)\n\ncolorbar = cm.ScalarMappable(cmap = cm.jet)\ncolorbar.set_array(color_map)\n'

In [301]:
def graficas (fig  , N=1, r = a, th=Th, phi=Ph, Rango=4, func=V, tag="$V$"):
    
    contador = 0. 
    for i in range(1,N+1):
        
        R = i*(r/N)
        progreso = (25./N)

        contador += progreso
        print(contador, "%")
        X = R*np.sin(Th)*np.cos(Ph)
        Y = R*np.sin(Th)*np.sin(Ph)
        Z = R*np.cos(Th)
        
       
        if type(func) == type(th):          
            color_map= np.real(func)
        else:
            color_map= np.real(func(R,th,phi,Rango))


        #color_map= np.real(func(R,th,phi,Rango))
        norm = colors.Normalize()
        color = norm(color_map)
        
        contador += progreso
        print(contador, "%")
        
        colorbar = cm.ScalarMappable(cmap = cm.jet)
        colorbar.set_array(color_map)
        

        if N != 2:
            ax = fig.add_subplot(np.ceil(N/2),np.ceil(N/2), i, projection="3d")
        else:
            ax = fig.add_subplot(N,1, i, projection="3d")
        
        contador+= progreso
        print(contador, "%")
        
        ax.set_xlim(-1,1)
        ax.set_ylim(-1,1)
        ax.set_zlim(-1,1)
        ax.set_title(r"Grafica de %s $(r, \theta, \phi)$ para $r = %2.2f$ "%(tag,R))
        ax.plot_surface(X, Y, Z,  rstride=1, cstride=1, facecolors=cm.jet(color))
        
        contador += progreso
        print(contador, "%")
        
        plt.colorbar(colorbar, shrink=0.7)
            
    print ("Puntos graficados: ",Th[0].size)
    print("Grado maximo del polinimio: ", Rango)
    print ("**************************************")

El campo eléctrico $\vec{E}$ se define como: $$ - \vec{\nabla} V = \vec{E}$$ Donde el operador gradiente debe incluirse en coordenadas esfericas, siendo: $$\vec{\nabla} = \frac{\partial}{\partial r} \hat{r} + \frac{1}{r} \frac{\partial}{\partial \theta} \hat{\theta} + \frac{1}{r~sin\theta} \frac{\partial}{\partial \phi} \hat{\phi}$$ Implementando lo anterior en una funcion, haciendo cortes para diferentes radios:


In [302]:
def E (R, TT, PP ,Rango=4):
    et,ep = np.gradient(np.real(V(R,TT,PP, Rango=Rango)))
    
    Et = (-1./(a*np.ones_like(TT)))*et                 #Factores de escala
    Ep = (-1./a*np.sin(TT))*ep
    return Et, Ep

Para hacer cortes en la coordenada radial, será necesario tomar $r =$ constante. Así, serán consideradas solo las componentes $E_{\theta}$ y $E_{\phi}$.


In [303]:
E_theta = E(a,Th,Ph)[0]
E_phi = E(a,Th,Ph)[1]

Graficando para $E_{\theta}$, para $r = a$:


In [307]:
figT= plt.figure(121)
graficas(figT, func=E_theta, Rango= 3, tag="$Et$")


25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  3
**************************************

Analogamente para $E_{\phi}$, para $r = a$:


In [309]:
figP= plt.figure(122)
graficas(figP, func=E_phi, Rango= 5, tag="$Ep$")


25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  5
**************************************

A continuación, una inocua funcion para guardar imagenes .png y hacer un gif despues.


In [311]:
def guardar (cuadros = 10, funcion=E_phi, archivo="Ep"):
    for i in range(1,cuadros+1):
        t = 1 - (i/float(cuadros))
        r = a - t 
        fig=plt.figure()
        graficas(fig, r = (a-t), func=funcion, tag="$E$")
        fig.savefig(archivo+"r=%1.2f_%d"%(r,(i-1))+".png")
        print(" --------->  Guardada %s"%(archivo+"_r=%d_%i"%(r,(i-1)) ))
        fig.clf
        plt.close

In [306]:
guardar()


25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_0
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_1
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_2
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_3
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_4
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_5
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_6
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_7
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=0_8
25.0 %
50.0 %
75.0 %
100.0 %
Puntos graficados:  100
Grado maximo del polinimio:  4
**************************************
guardada Ep_r=1_9

In [306]: