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]:
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]:
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$")
Analogamente para $E_{\phi}$, para $r = a$:
In [309]:
figP= plt.figure(122)
graficas(figP, func=E_phi, Rango= 5, tag="$Ep$")
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()
In [306]: