Función de Green para el problema de Dirichlet interior a dos esferas concéntricas conectadas a tierra.

Los parámetros son "a" para el radio menor, "b" para el radio mayor y la carga puntual unitaria está ubicada en r=d

Par = {a -> 2, b -> 4, 
          d -> 3.5};
x = (b/a);
A = (a/d);
B = (b/d);

Escribo la función de Green tomando 30 términos del desarrollo en serie...

f[r_, α_] := 
        Sum[((1/x)^i ((1/(Sqrt[(r^2 + ((d/x)^(2  i))^2 - ((2 r d)/x)^(2  i) Cos[α])])) - 
                          (A/(Sqrt[(r^2 + (((A^2 d)/x)^(2  i))^2 - ((2 r d A^2)/x)^(2  i) Cos[α])]))) + 
                    x^i ((1/(Sqrt[(r^2 + (x^(2  i) d)^2 - 2 r d x^(2  i) Cos[α])])) - 
                          B/(√(r^2 + (B^2 x^(2  i) d)^2 - 
                                    2 r d x^(2  
                                        i) B^2 Cos[α])))) /. 
                Par, {i, 0, 30}] - 
            (1/√)(r^2 +  d^2 - 2 r d Cos[α]) /. 
          Par;

Grafico la función barriendo ángulos de la forma π/10^n . A medida que n tiende a infinito nos acercamos a la posición angular de la carga puntual (observar el comportamiento en la vecindad de r=3.5 que es la ubicación de la carga real...)

Table[
      Plot[f[r, π/(10^n)], {r, 2, 4}, PlotRange -> All], {n, 0, 
        4}]

Ahora graficamos la función de Green en una lista de gráficos tridimensionales para ver su variación característica a distintas escalas. Observar que la condición de contorno solo se pone en evidencia a escalas pequeñas mientras que la presencia de la carga real aparece evidente cuando aumentamos la escala....

CartesianForm[u_, v_] = 
        f[r, α] /. {r -> (Sqrt[u^2 + v^2]), 
            Cos[α] -> 
              (u/(Sqrt[u^2 + v^2]))};
Green[u_, v_] := 
        If[(a /. Par) < (Sqrt[u^2 + v^2]) < (b /. Par), 
          CartesianForm[u, v], 0];
Table[
        Plot3D[Green[u, v], {u, -5, 5}, {v, -5, 5}, 
          PlotPoints -> 80, 
          PlotRange -> {0, 10^(-n)}], {n, -1, 
          6}];

Ahora vemos un gráfico de contornos de la función de Green en escala logarítmica para eliminar las diferencias de variaciones características según la escala.

Green2[u_, v_] := 
        If[(a /. Par) < (Sqrt[u^2 + v^2]) < (b /. Par), 
          Log[CartesianForm[u, v]], 0];
ContourPlot[Green2[u, v], {u, -5, 5}, {v, -5, 5}, 
        PlotPoints -> 80, ContourLines -> False];

Carga inducida en las superficies de las esferas:

Sigmab y Sigmaa son las cargas inducidas en las superficies de las esferas (observar que, como deben, ambas densidades son negativas)

Sigmab = (1/(4  π)) 
            D[f[r, α], 
              r] /. {r -> (b /. Par)};
Plot[Sigmab, {α, -π, π}, 
      PlotRange -> All]

Sigmaa = -(1/(4  π)) 
            D[f[r, α], 
              r] /. {r -> (a /. Par)};
Plot[Sigmaa, {α, -π, π}, 
      PlotRange -> All]