function returning 0

Question

Why is my function returning 0?

In[95]:= n = 25;
r[x_, y_] = Sqrt[x^2 + y^2];
\[Theta][x_, y_] = ArcTan[x, y];
f[y] = Piecewise[{{100, 0 <= y <= 1}, {0, -1 <= y <= 0}}];

In[109]:= 
c = Table[(2*l + 1)/2*
    Integrate[f[y]*LegendreP[l, 0, y], {y, -1, 1}], {l, 0, n}];

In[100]:= 
T[x_, y_] = 
 Parallelize[
  Sum[c[[n]]*r[x, y]^l*LegendreP[l, 0, Cos[\[Theta][x, y]]], {l, 0, 
    n}]]

Out[100]= 0

I don't understand why T[x,y] is returning only 0.

---

Edit 2:

So if I start l from 1 and add in c[[0]] manual to the sum the code works but something is still going wrong.

The top of northern hemisphere is supposed to be held at 100 and the southern hemisphere is supposed to be held at 0.

When I plot this for different n, the temperature distribution switches to 100 on the top to 100 on the bottom and similarly for 0.

And the 100 on the top is always chopped out.

This is n = 25 and n = 23, respectively.

DensityPlot[T[x, y], {x, -1, 1}, {y, -1, 1}, 
 ColorFunction -> "Rainbow", 
 RegionFunction -> Function[{x, y}, 0 < x^2 + y^2 < 1]]

What am I doing wrong? Why is 100 chopped off, why is it flipping position at different n, and why do I have to add in c[[0]] instead of letting it integrate l = 0 and summing at l = 0?

If n is too large, then the image is:

Edit 3:

Here is the code I am using. Try n = 25, n = 35, and n = 50 to see what happens.

n = 35;
r[x_, y_] := Sqrt[x^2 + y^2];
\[Theta][x_, y_] := ArcTan[y, x];
f[x_] := Piecewise[{{100, 0 <= x <= 1}, {0, -1 <= x <= 0}}];

c = Table[(2*l + 1)/2*
    Integrate[f[x]*LegendreP[l, 0, x], {x, -1, 1}], {l, 1, n}];

T[x_, y_] = 
  50 + Sum[c[[n]]*r[x, y]^l*LegendreP[l, 0, Cos[\[Theta][x, y]]], {l, 
     0, n}];

Plot3D[T[x, y], {x, -1, 1}, {y, -1, 1}, Mesh -> None, Boxed -> False, 
 ColorFunction -> "Rainbow", 
 RegionFunction -> Function[{x, y}, 0 < x^2 + y^2 < 1], 
 PlotRange -> All]

DensityPlot[T[x, y], {x, -1, 1}, {y, -1, 1}, 
 ColorFunction -> "Rainbow", 
 RegionFunction -> Function[{x, y}, 0 < x^2 + y^2 < 1], 
 PlotRange -> All]

Answer 1

Here's what I think you probably want:

n = 50;
r[x_, y_] := Sqrt[x^2 + y^2];
θ[x_, y_] := ArcTan[y, x];
f[x_] := Piecewise[{{100, 0 <= x <= 1}, {0, -1 <= x <= 0}}];

c = Table[(2 * l + 1) / 2 * Integrate[f[x]*LegendreP[l, 0, x], {x, -1, 1}], {l, 0, n}];

T[x_, y_] = Sum[c[[1 + l]] * r[x, y]^l * LegendreP[l, 0, Cos[θ[x, y]]], {l, 0, n}];

Then:

Plot3D[T[x, y], {x, -1, 1}, {y, -1, 1}, Mesh -> None, Boxed -> False, 
 ColorFunction -> "Rainbow", 
 RegionFunction -> Function[{x, y}, 0 < x^2 + y^2 < 1], 
 PlotRange -> All]

DensityPlot[T[x, y], {x, -1, 1}, {y, -1, 1}, 
 ColorFunction -> "Rainbow", 
 RegionFunction -> Function[{x, y}, 0 < x^2 + y^2 < 1], 
 PlotRange -> All]

For the sake of speed, consider using Dot and the listability of Power and LegendreP:

T2[x_?NumericQ, y_?NumericQ] := 
  With[{r0 = r[x, y]}, c.(r0^Range[0, n]*LegendreP[Range[0, n], 0, y/r0])];

It's more than twice as fast:

DensityPlot[T[x, y], {x, -1, 1}, {y, -1, 1}, 
  ColorFunction -> "Rainbow", 
  RegionFunction -> Function[{x, y}, 0 < x^2 + y^2 < 1], 
  PlotRange -> All] // AbsoluteTiming
DensityPlot[T2[x, y], {x, -1, 1}, {y, -1, 1}, 
  ColorFunction -> "Rainbow", 
  RegionFunction -> Function[{x, y}, 0 < x^2 + y^2 < 1], 
  PlotRange -> All] // AbsoluteTiming

One can also compute a formula for the coefficients (it takes about 0.5 sec.):

coeff[l_] = (2*l + 1) / 2 * Integrate[f[x] * LegendreP[l, 0, x], {x, -1, 1}]
(* (25 (1 + 2 l) Sqrt[π])/(Gamma[1 - l/2] Gamma[(3 + l)/2]) *)

Then computing c can be quite fast:

n = 100;
c = Table[(2*l + 1)/ 2 * Integrate[f[x] * LegendreP[l, 0, x], {x, -1, 1}], {l, 0, n}]; //
      AbsoluteTiming
c = coeff /@ Range[0, n]; // AbsoluteTiming

(* {8.634086, Null} *)
(* {0.002338, Null} *)