Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]]

enter image description here

enter image description here

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:

enter image description here

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]
share|improve this question
1  
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. –  b.gatessucks Oct 5 '13 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 :. –  Sektor Oct 5 '13 at 7:49
    
good practice, but it should work since he only uses f with a literal y argument. –  george2079 Oct 5 '13 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.. –  george2079 Oct 5 '13 at 12:26
1  
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'. –  Ymareth Oct 5 '13 at 17:40

1 Answer 1

up vote 3 down vote accepted

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]

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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} *)
share|improve this answer
    
Yes, that is what I was striving for. Is it possible to plot a 3D sphere with this correct gradient as well? –  dustin Oct 5 '13 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. –  Michael E2 Oct 5 '13 at 21:03
    
I mean since this a solution to Laplace on sphere can the sphere with the heat map be plotted? –  dustin Oct 5 '13 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. –  dustin Oct 5 '13 at 21:10
1  
@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 ]' –  george2079 Oct 7 '13 at 17:35

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.