0
$\begingroup$

The problem I'm having is that the function definition given below continually gives me baffling error messages.

I'm trying to turn a list of 3D points into a function that approximates a curve by connecting consecutive points with straight lines. If anyone has a smarter way to do that than what I describe below (a smooth curve would be fine too), I'll accept suggestions! Note that it's vital for the curve to be parametric.

The definition says, in words: "define each location on the 'curve' as a weighted average of the two nearest given points, where the points are provided in a list, and each element in the list is a list of three coordinates. (e.g. {{1,1,1},...})." The effect of weighting (given by 'frac' and '1-frac') is to effectively construct a continuous straight line between each two given points. The 'If' statement ensures that when constructing the last point, the function won't attempt to look for a point after the last point.

Here's the definition:

curve[s_] := 
  (frac = 1 - Mod[s*100, 1]) points[[low = Floor[s*100] + 1]] + 
    If[s < 1.0, (1 - frac) points[[low + 1]], 0]

curve3D = ParametricPlot3D[curve[s], {s, 0.0, 1.0}]

The call to ParametricPlot3D gives me this error message:

Part::pspec: Part specification 1+Floor[100 s] is neither a machine-sized integer nor a list of machine-sized integers. >>

or sometimes (after a few minutes of debugging) it starts giving this instead:

Part::partw: Part 102 of {{0.,2.0944,15.},{0.015708,2.09432,14.9963},{0.0314159,2.09411,14.9852},{0.0471239,2.09375,14.9667},<<44>>,{0.753982,1.94388,7.97093},{0.76969,1.93806,7.73558},<<51>>} does not exist. >>

The first set of messages seem to be saying that 's' isn't getting evaluated, which makes no sense to me.

The second set say that I'm trying to access the 102nd point (there are only 101 points), which I believe I'm avoiding with the If statement. The message even persists if I change the parametric range of s to {s, 0.0, 0.5} and if I clear the variable "low" beforehand.

Can anyone explain either of these error messages, and what I'm doing wrong? Is it a misuse of ParametricPlot3D?

For further clarity, here's the output of curve3D (rotated on its side to save space). Note the odd gap in the middle. As you change the parametric range of s (e.g. {s, 0.0, 0.8}) the number and location of the gaps changes. Super weird stuff.

graph

Also, if you want to run my code yourself, here's how you construct the list 'points':

potential[x_, y_] := 20 (Cos[x]^2 + Cos[y]^2 + Cos[x] Cos[y])
nmax = 100;
points = {};

For[n = 0, n <= nmax, n++,
  xn = N[(π/2)*n/nmax];
  ptrule = FindMinimum[potential[xn, y], {y, π/2}];
  pt = {xn, y, potential[xn, y]} /. ptrule[[2]];
  AppendTo[points, pt];]
$\endgroup$
  • $\begingroup$ I can't seem to reproduce your error messages when running the first portion of your code (the ParametricPlot3D). Have you tried restarting the kernel on the off chance that some stray definitions are giving you trouble? $\endgroup$ – MarcoB Jul 10 '15 at 19:11
  • $\begingroup$ Yeah I've restarted the kernel many a time now. :/ $\endgroup$ – maxdp Jul 10 '15 at 19:24
  • 1
    $\begingroup$ I'm afraid that I have to repeat the same observation: I copied and run your code into a fresh notebook, and received no errors. I obtain more or less the same output from ParametricPlot3D you showed. I am using MMA v. 10.1.0 on Win7-64bit. What version are you on, just in case? I just can't seem to reproduce your problem. Can you simplify your code further to a minimal example that still produces that error? $\endgroup$ – MarcoB Jul 10 '15 at 23:53
  • $\begingroup$ @MarcoB I confirm that the code runs fine using MMA v. 10.1.0 on Win8.1-64bit. $\endgroup$ – bbgodfrey Jul 11 '15 at 17:13
  • 1
    $\begingroup$ The code runs without error messages using MMA v 10.1 on Mac OS 10.10.4; however, the gap appears. The gap can be removed by setting the PlotPoints option for ParametricPlot3D to a large value. $\endgroup$ – Bob Hanlon Jul 11 '15 at 21:40
1
$\begingroup$

I've run this code on V8.04, V9.0.1, V10.0.2 and V10.1 (Windows 8.1) and I see the first error reported by the OP in v9.0.1 only.

In V9.0.1 you can prevent the error from occurring by changing the definition of curve from curve[s_] to curve[s_?NumericQ]. A related answer here.

I do not see the second error message in any of the four versions.

The gap can be removed by setting the option PlotPoints to for instance 1000:

Mathematica graphics

Interestingly, in V9 after using curve[s_?NumericQ] there's no gap.

$\endgroup$
  • $\begingroup$ Thank you, this was it. I'm working in V9.0.0. I've never heard of ?NumericQ so thanks for pointing me to this (bizarre) fix. Guess I'll be upgrading my version soon... $\endgroup$ – maxdp Jul 13 '15 at 16:23
  • 1
    $\begingroup$ And thank you for the accept. ?NumericQ is less bizarre than you think. In many functions (like minimization functions) Mathematica first does a symbolic call to the function it is given as an argument. This goes wrong if the (user) function does not expect symbolic arguments. Using ?NumericQ to restrict the input to numeric types prevents this. $\endgroup$ – Sjoerd C. de Vries Jul 13 '15 at 20:30
2
$\begingroup$

If curve is given a non-numeric input the Part specification is not valid and we get:

curve["atom"]

Part::pkspec1: The expression 1+Floor[100 atom] cannot be used as a part specification. >>

Since symbolic analysis may be used by many Mathematica functions, including plotting, you should restrict the input of your curve function to numeric values using pattern matching.

The second error may be reproduced by entering a s value that is slightly too high:

curve[1.01]

Part::partw: Part 102 of {{0.,2.0944,15.},{0.015708,2.09432,14.9963},<<47>>,{0.76969,1.93806,7.73558},<<51>>} does not exist. >>

Plot functions may try values slightly outside the given plot range. This can cause problems just such as this, but I seem to recall that there is a reasonable explanation for this design choice; if I find a reference I shall link it. A solution is to further restrict the input of curve.

Also you should be localizing frac and low with Module.

Combined:

Clear[curve]

curve[s_] /; 0 <= s <= 1 :=
  Module[{frac, low},
    (frac = 1 - Mod[s*100, 1]) points[[low = Floor[s*100] + 1]] + 
      If[s < 1.0, (1 - frac) points[[low + 1]], 0]
  ]
$\endgroup$
1
$\begingroup$

This answer is not concerned with the error messages your code produces, but to suggest your problem could have been avoided by doing more up-front analysis with Mathematica.

I will begin by looking at the behavior of your potential function.

potential[x_, y_] := 20 (Cos[x]^2 + Cos[y]^2 + Cos[x] Cos[y])
Plot3D[potential[x, y], {x, 0, π/2}, {y, 0, π/2}, AxesLabel -> Automatic]

plot3d

You are seeking the curve generated by the points {x, y, potential[x, y]}, where y minimizes the single variable function potential[x0, y] for some 0 <= x0 <= π/2. The above plot suggests that this curve is not a valley running through the 3D potential surface, but the intersection of the surface with the bounding box wall at x = π/2. Another plot will make this even clearer. It is a 2D plot of the curves of potential for fixed x values. That is, we are looking at the potential surface's intersection with planes parallel to the xz-plane at various x distances.

Plot[Evaluate[Through[#2[y]]], {y, 0, π/2}, 
  PlotLegends -> LineLegend[Row[{"x = ", #}] & /@ #1]] & @@ 
  Transpose @ Table[{x, With[{x = x}, potential[x, #] &]}, {x, 0, π/2, π/8}]

slices

Although I am making a heuristic argument, I find it convincing, and I assert that a plot of the curve you are looking for is given by the following simple expression:

    ParametricPlot3D[{x, π/2, potential[x, π/2]}, {x, 0, π/2}]

curve

$\endgroup$
  • $\begingroup$ Thanks for the feedback, but what I'm getting is actually a valley. It'd be easier to see if you used the range {x, 0, pi}, {y, 0, pi}, where the entire valley lies. $\endgroup$ – maxdp Jul 13 '15 at 16:32
  • $\begingroup$ @maxdp. I used the range you gave in your question. If you were actually using a different range than your question is ill-posed. $\endgroup$ – m_goldberg Jul 13 '15 at 18:12
  • $\begingroup$ I don't see where in my question I mentioned that range. The only place the range is given is in my diagram, which has axes showing {x, 0, pi/2}, {y, 2ish, pi/2}, {z, 15ish, 0}. And in FindMinimum at the bottom I guess that y is near pi/2, which it is. Maybe I should've mentioned the range more explicitly, but it didn't seem important (and ultimately it wasn't). $\endgroup$ – maxdp Jul 13 '15 at 18:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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