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.

I've been having some trouble with Plot3D. I have a quite complex function of two variables (actually three but the first one doesn't matter) that returns a list with two elements. Then I wrote the following code to see the behaviour of the first element:

Module[{R = 1},
 Plot3D[
  g[R, θ, z][[1]],
  {θ, 0, Pi}, {z, 0, 4}, 
  AxesLabel -> {θ, z, Subscript[F, θ]}
  ]
 ]

Well the thing is: when the coordinate $\theta=0$, the first element is zero and the plot doesn't show that:

enter image description here

I know that the problem isn't from the function because if I plot the first element keeping $\theta=0$:

Module[{R = 1},
 Plot3D[
  g[R, 0, z][[1]],
  {θ, 0, Pi}, {z, 0, 4}, 
  AxesLabel -> {θ, z, Subscript[F, θ]}
  ]
 ]

(in the code I only changed that $\theta \to 0$) it outputs correctly:

enter image description here

Moreover, If I use Manipulate it behaves as expected; so I must be doing something wrong with Plot3D... can anyone help me?

PS: Here is the horrible definition of the function $g$:

g[R_, \[Theta]_, 
   z_] := (Module[{For1 = 0, For2 = 0, d1 = 0, d2 = 0, Forca = {0, 0},
      nmax = 10},

    If[\[Theta] == 0,
     For[n = 0, n < nmax,
       d1 = Sqrt[R^2*(2 Pi*n)^2 + z^2];
       If[n == 0,
        If[d1 + 0.0 < R*Pi,
          For1 = -1/(2 Pi*d1^2);
          Vec1 = {0, z};
          ,
          For1 = -1/(4*d1^2*ArcSin[R*Pi/d1]);
          Vec1 = {0, z};
           ];
        ,
        For1 = -1/(4*d1^2*ArcSin[R*Pi/d1]);
        Vec1 = {0, z};
        For2 = -1/(4*d1^2*ArcSin[R*Pi/d1]);
        Vec2 = {0, z};
         ];
       Forca += For1*Vec1 + For2*Vec2;
       ++n;
        ];
       (*Here works iff \[Theta]=0*)
       ,
       For[n = 0, n < nmax,
         d1 = Sqrt[R^2*(\[Theta] + 2 Pi*n)^2 + z^2];
         d2 = Sqrt[R^2*(-2 Pi*(n + 1) + \[Theta])^2 + z^2];

         If[d1 + 0.0 < R*Pi,
             For1 = -1/(2 Pi*d1^2);
             Vec1 = {\[Theta], z};
             ,
              For1 = -1/(4*d1^2*ArcSin[R*Pi/d1]);
              Vec1 = {\[Theta] + 2 Pi*n, z};
              ];

          If[d2 + 0.0 < R*Pi,
              For2 = -1/(2 Pi*d2^2);
              Vec2 = {-2 Pi*(n + 1) + \[Theta], z};
              ,
              For2 = -1/(4*d2^2*ArcSin[R*Pi/d2]);
              Vec2 = {-2 Pi*(n + 1) + \[Theta], z};
               ];
           Forca += For1*Vec1 + For2*Vec2;
           ++n;
           ]
     (*Here works for all points except \[Theta]=0*)
       ];
         N[Forca]
    ]
   );
share|improve this question
    
@Jens I was answering a question in another comment, now deleted. Deleting –  belisarius Mar 10 '13 at 19:40
add comment

1 Answer

up vote 4 down vote accepted

This is probably due to the degenerate nature of your function.

Take for example this slice (z = 1):

Plot[g[1, θ, 1][[1]], {θ, -2, 2}, PlotRange -> All]

Mathematica graphics

It looks like it's continuous a θ=0 but it isn't:

g[1, #, 0.5][[1]] & /@ {-0.0001, 0, -0.0001}

{0.07960574141, 0., 0.07960574141}

Things gets worse for lower z values:

Plot[g[1, θ, 0.01][[1]], {θ, -2, 2}, PlotRange -> All]

Mathematica graphics

The zero value at θ=0 is probably hard to find with most samplings of the parameter space.

A workaround would be to use ListPlot3D:

Module[{R = 1}, 
 ListPlot3D[
  Flatten[Table[{θ, z, g[R, θ, z][[1]]}, {z, 0.001, 4, 0.1}, {θ, 0, Pi, Pi/40}], 1], 
  AxesLabel -> {θ, z, Subscript[F, θ]}]
]

Mathematica graphics

share|improve this answer
add comment

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.