Quit[]

(*Paquete para correr*)

Needs["PlotLegends`"]

(* Consideremos las soluciones que encontramos para el problema que hicimos en clase. El problema consistía en un plano ubicado en z=0 que estaba completamente a tierra excepto por un circulo de radio "a" centrado que estaba a potencial V0. Con el metodo de Green pudimos obtener el potencial en el semi-eje z>0 y en el semi-eje z<0. Luego, usando la extensión analítica pudimos obtener el potencial en todo el espacio como una serie en "r", con polinomios de Legendre del Cos[θ]. Por otra parte vimos también que el mismo problema podía ser atacado usando separación de variables en coordenadas cilíndricas obteniendo una expresión del mismo como una transformada de Hankel *)

Definimos el potencial que obtuvimos por el primer método. Este está dado por una serie infinita que decidimos cortar en el valor “n”. A su vez, como tenemos simetría de revolución, nos interesará graficar el potencial en una superficie fija (e imaginar que su forma tiene simetría de revolución) y elegimos el plano x=0. Así, r=(Sqrt[(y^2)+(z^2)])y Cos[θ]=(z/(Sqrt[(y^2)+(z^2)])). Por último, elegimos fijar V0=1 y a=1 de manera que medimos el potencial en unidades de V0 y medimos las distancias en unidades de “a”. Como vimos en la práctica debemos definir el potencial por tramos separando para r<a y r>a :

ϕ1[y_,z_,n_]:=1-1Sign[z]Sum[((Pochhammer[1/2,l](((-1))^l))/(l!))((((Sqrt[(y^2)+(z^2)])))^(2l+1))LegendreP[2l+1,(z/(Sqrt[(y^2)+(z^2)]))],{l,0,n}]/;(y^2)+(z^2)<1;
ϕ1[y_,z_,n_]:=-1Sign[z]Sum[((Pochhammer[1/2,l](((-1))^l))/(l!))((((1/(Sqrt[(y^2)+(z^2)]))))^(2l))LegendreP[2l-1,(z/(Sqrt[(y^2)+(z^2)]))],{l,1,n}]/;(y^2)+(z^2)≥1

Podemos entonces graficar la intensidad del potencial en el plano x=0. Vemos que en z=0 el potencial es nulo para r>a y muy intenso para r<a. Visto desde lejos r>>a el potencial es parecido al “doble hongo” de un doble dipolo (como veremos más abajo).

ShowLegend[DensityPlot[ϕ1[y,z,10],{y,-3,3},{z,-3,3},ColorFunction->"SunsetColors"],{ColorData["SunsetColors"][1-#1]&,10," Max","Min",LegendPosition->{1.1,-.4}}]

Si graficamos la intensidad del potencial cerca del plano z=0 y en bien dentro de la región interior (r<a) obtenemos el potencial de un plano infinito:

ShowLegend[DensityPlot[ϕ1[y,z,10],{y,-1/4,1/4},{z,-1/4,1/4},ColorFunction->"SunsetColors"],{ColorData["SunsetColors"][1-#1]&,10," Max","Min",LegendPosition->{1.1,-.4}}]

Aquí definimos el modo k-esimo de la solución exterior r>a (y seteamos a cero la solución interior r<a porque solo nos interesa el desarrollo multipolar exterior). Como vemos de arriba, las contribuciones son todas impares en k. Luego graficamos con “Manipulate” para ver las contribuciones de cada modo. Es interesante ver el patron de nodos y máximos que surge al aumentar k. Esta claro de este análisis, al comparar con el gráfico de arriba, que es el modo k=1 el que domina la forma de la solución exterior.

ϕmode[y_,z_,k_]:=0/;(y^2)+(z^2)<1;ϕmode[y_,z_,k_]:=-1Sign[z]((Pochhammer[1/2,((k+1)/2)](((-1))^((k+1)/2)))/((((k+1)/2))!))((((1/(Sqrt[(y^2)+(z^2)]))))^(k+1))LegendreP[k,(z/(Sqrt[(y^2)+(z^2)]))]/;(y^2)+(z^2)≥1

Manipulate[ShowLegend[DensityPlot[ϕmode[y,z,k],{y,-2,2},{z,-2,2},PlotPoints->100,ColorFunction->"SunsetColors"],{ColorData["SunsetColors"][1-#1]&,10," Max","Min",LegendPosition->{1.1,-.4}}],{k,1,11,2}]

Aquí en cambio definimos nuestro potencial como una integral númerica que aproxima al potencial que obtuvimos como una transformada de Hankel. Vemos que este gráfico es computacionalmente mucho más costoso que la serie de polinomios de Legendre.

ϕ2[y_,z_]:=NIntegrate[BesselJ[1,k]BesselJ[0, k Abs[y]](E^(-k Abs[z])),{k,0,∞}]

p2=ShowLegend[DensityPlot[ϕ2[y,z],{y,-3,3},{z,-3,3},ColorFunction->"SunsetColors"],{ColorData["SunsetColors"][1-#1]&,10," Max","Min",LegendPosition->{1.1,-.4}}]

Aquí graficamos ambos resultados a la par. Esto es, el resultado como una serie en “r” con polinomios de Legendre y el resultado como una transformada de Henkel. Vemos la similitud entre ambos gráficos.

GraphicsRow[{DensityPlot[ϕ1[y,z,10],{y,-3,3},{z,-3,3},ColorFunction->"SunsetColors"],DensityPlot[ϕ2[y,z],{y,-3,3},{z,-3,3},ColorFunction->"SunsetColors"]}]

El gráfico izquierdo a continuación es el doble hongo que me refería más arriba. Corresponde al campo de un dipolo apuntando hacia arriba (si lo miro desde arriba) y un dipolo apuntando hacia abajo (si lo miro desde abajo). Grafico a la derecha la funcion ϕ1[y,z] (nuestra primer solucion) con una escala que va -10 a 10 en ambas variables, como si estuviera haciendo un “zoom back” de nuestra solución. Vemos que ambos gráficos son muy similares.

GraphicsRow[{DensityPlot[((Abs[z])/((((Sqrt[(y^2)+(z^2)])))^3)),{y,-10,10},{z,-10,10},PlotPoints->100,ColorFunction->"SunsetColors"],DensityPlot[ϕ1[y,z,10],{y,-10,10},{z,-10,10},ColorFunction->"SunsetColors",PlotPoints->100]}]