# function returning 0

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]

• Regarding the first part, you use c[[n]] which happens to be zero; maybe you meant to use c[[l]] in the summation. Also you need to define f[y_] with an underscore. Commented Oct 5, 2013 at 7:45
• Adding to b's comment you should define your functions like this f[x_] := r.h.s. As you can see you have omitted the :. Commented Oct 5, 2013 at 7:49
• good practice, but it should work since he only uses f with a literal y argument. Commented Oct 5, 2013 at 12:18
• i dont know if it pertains to the problem, but this seems a poor use of parallelize -- the integration is already done when you define c. If anything you should use paralletable on the c definition.. Commented Oct 5, 2013 at 12:26
• f[y]='something involving y' will not work (as expected) in Table[something involving f[y],{y,some range}] as the y inside table is not the global symbol y. Try f[y] = 4 and Table[f[y], {y, 1, 3}] and you get {f[1], f[2], f[3]} not {4, 4, 4}. You do need f[y_]= 'something involving y' or more likely f[y_]:= 'something involving y'. Commented Oct 5, 2013 at 17:40

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} *)

• Yes, that is what I was striving for. Is it possible to plot a 3D sphere with this correct gradient as well? Commented Oct 5, 2013 at 20:46
• @dustin Forgive me, but I can't figure out what gradient you're referring to or how it would map onto a sphere. But I'm pretty sure that if you know the mathematical formulation, then it will be possible to do. Commented Oct 5, 2013 at 21:03
• I mean since this a solution to Laplace on sphere can the sphere with the heat map be plotted? Commented Oct 5, 2013 at 21:05
• Also, is there a way to speed up the finding of c? When I increase n, the calculation of c slows down immensely. Commented Oct 5, 2013 at 21:10
• @dustin ParallelTable speed it up quit a lot (Michaels approach wont always work for a general not-piecewise-constant f..). If you retain only floating values that helps as well.. 'ParallelTable[ ... Integrate[] //N ]' Commented Oct 7, 2013 at 17:35