Mathematica - sorting NSolve solutions

[Dale]

At the points where it jumps, why is it jumping? Also, is there any rule about the regions where the roots come from, i.e. I know that the order changes but can you say that there is always exactly one root with y>.2, exactly one with -.2<y<0 and 0<x<1, ...

 

[Bill Simpson]

If you post the 5-6 data files and your latest notebook that did the six color plot above along with some directions on how to use your notebook and what you really want as output then I'll look at it and see if I can get something to work.

 

[MartinV]

At the points where it jumps, why is it jumping?

If I'd known that, I would be done with my bachelor's degree by now. It's got something to do with the fact that the lower, apple-shapped formation, is moving in the opposite direction as the upper one (hat-like).

The inputs are traveling from left to right. The 'hat formation' series (cyan first, then magenta, then black, magenta again, cyan in the end) also goes from left to right. The lower series starts at (1,0) (magenta), then shifts to blue and cyan. At the bottom another picture forms (green), cyan elipses disappear, then the green turns into blue and black.
The top of the "apple" is a separate picture.

Bill Simpson, I usually run the code at the beginning of this thread with input files. I can zip those files and send them to you. Then all you need to do is run the program (it runs approx 15 minutes with 187 input files.)

 

[Dale]

In post 13 you gave 10 inputs. I was able to use your code to replicate one group of 5 ellipses, all neatly separated. As an intermediate to zipping up the whole input, perhaps you could find one spot where it "jumps" and give the 10 inputs before and the 10 inputs after. Then maybe I could figure out why it is jumping. You could do that either with the original jumps without my sorting function or with the new jumps with my sorting function. Either way is fine, just let me know.

 

[MartinV]

I have 41 files ready. How do I send them?

 

[MartinV]

Ok, I managed to separate the ellipses using FindClusters[]. The code now looks like this:

SetDirectory["C:\..."];

p = 1.0; q = 0; eps = 0.5; M = 101;

For[j = 1, j <= M, j += 1,
P = Import["P" <> ToString[j] <> ".txt", "Table"]; b = {}; c = {}; d = {};
For[i = 1, i <= 200, i++,
a = NSolve[{x - (1 - eps)*x/(x^2 + y^2) - eps*(x - p)/((x - p)^2 + (y - q)^2) == P[[i, 1]],
y - (1 - eps)*y/(x^2 + y^2) - eps*(y - q)/((x - p)^2 + (y - q)^2) == P[[i, 2]]}, {x, y}]; b = Append[b, a];
] ; c = {x, y} /. Partition[Flatten, 2]; d = Select[c, FreeQ[#, Complex] &];
{e[j],f[j], g[j]} = FindClusters[d, 3, Method -> "Agglomerate"];
]

Now I have one more thing left to do. I have to calculate the surface area of those ellipses.

This is what I have so far:

surf[x_] :=
Do[k =
x[[i, 1]]*(x[[i + 1, 2]] - x[[i, 2]]) -
x[[i, 2]]*(x[[i + 1, 1]] - x[[i, 2]]), {i, 1, Length[x]}]

I was planning to use surf[] on every array I get out of FindClusters[] like this:

surf[e[55]] (55 is the number of a random source. Input file P55.txt gives me solutions e[55],f[55],g[55])

Somehow it's not working. Index i begins at 200 instead of 1, probably because I used i as index before. I thought by doing {i,1,Length[x]} it would reset it.

Anway, I'm trying to do the integral over those points. I need to calculate a cross product between coordinates of points and the distance between them. If anyone knows a faster way to do this, let me know.

 

[MartinV]

This is what I have now:

surf[x_] :=
k = Table[x[[i, 1]]*(x[[i + 1, 2]] - x[[i, 2]]) - x[[i, 2]]*(x[[i + 1, 1]] - x[[i, 2]]), {i, Length[x]}]

It still doesn't work and I don't know why. I really need help here.

 

[Bill Simpson]

It still doesn't work and I don't know why. I really need help here.

Restart the kernel.
Evaluate the whole notebook in order up to but just before where you are going to evaluate surf[].
Create two lines:

mydata=WhatEverItIs;
surf[mydata]

Delete everything else, evaluate those two lines and past the result here.
You want to guarantee that I can evaluate this and get exactly the same thing you get, without needing hundreds of files or thousands of data points or perhaps hundreds of lines of Mathematica that nobody but you has or knows what they do. Simplify it down to those two lines (actually WhatEverItIs just needs to be big enough to allow someone to discover what the problem really is, no more and no less).

I really need to see exactly what the details of your WhatEverItIs so I can try to tell you why it isn't working and what to try to do to hopefully fix it.

 

[MartinV]

My WhatEverItIs (the data I need to feed into surf) is 200 pairs of numbers. Part of it:

{{-1.00704, 0.433162}, {-1.00634, 0.436155}, {-1.00568,
0.43911}, {-1.00506, 0.442026}, {-1.00448, 0.444898}, {-1.00395,
0.447725}, {-1.00346, 0.450502}, {-1.00301, 0.453228}, {-1.00261,
0.455899}, {-1.00226, 0.458513}, {-1.00195, 0.461068}, {-1.00169,
0.463561}, {-1.00149, 0.465989}, {-1.00133, 0.468352}, {-1.00122,
0.470646}, {-1.00116, 0.472869}, {-1.00115, 0.47502}, {-1.00119,
0.477097}, {-1.00128, 0.479098}, {-1.00142, 0.481022}, {-1.00161,
0.482866}, {-1.00185, 0.484631}, {-1.00215, 0.486314}}

The output of surf[mydata] is:

0.496875, 0.501302, 0.50583, 0.510457, 0.515175, 0.519982, 0.524871, \
0.529837, 0.534875, 0.53998, 0.545145, 0.550365, 0.555634, 0.560946, \
0.566296, 0.571677, 0.577084, 0.58251, 0.587949, 0.593396, 0.598843, \
0.604286, 0.609718, -0.430136 (-0.430136 + {{-1.00704,
0.433162}, {-1.00634, 0.436155}, {-1.00568, 0.43911}, {-1.00506,
0.442026}, {-1.00448, 0.444898}, {-1.00395, 0.447725}, {-1.00346,
0.450502}, {-1.00301, 0.453228}, {-1.00261,
0.455899}, {-1.00226, 0.458513},

then in the middle I get this:

{-1.00779, 0.430136}}[[201, 1]]) - 1.00779 (-0.430136 + {{-1.00704, 0.433162}, {-1.00634, 0.436155},

Those intial numbers are probably what I'm looking for but I want a plain array of numbers. 200 lines in input, 200 lines in output. In truth I need the sum of those 200 output numbers but for now I would be happy if I just get those 200.

 

[Bill Simpson]

Ah, good, we are finally making some progress.

You have a vector of 200 {x,y}.
Your surf has Table[x[[i, 1]]*(x[[i + 1, 2]] - x[[i, 2]]) - x[[i, 2]]*(x[[i + 1, 1]] - x[[i, 2]]), {i, Length[x]}].

Notice that i ranges over the length of your vector BUT your Table has x[[i+1,2]] and x[[i+1,1]] in it. Those are going to "subscript off the end" of your vector. Unless you have turned it off (you didn't actually show everything I asked for, input, output, error messages, etc, all the real raw stuff to help diagnose your problem) you should have gotten some error, perhaps like
Part::partw: "Part 201 of {YourHugeXList} does not exist.

So I don't know if this is really going to fix it, but you can have {i,Length[x]-1} so your Table doesn't try to reach beyond the end of your list. Or maybe you need to carefully think about your subscripting so you really use all the points and nothing beyond either end.

Does this help? See if you can fix your surf[] from this and if that doesn't fix everything then post another message with all the input and output and errors and what is wrong so we can track this down. And manually check that output data to make sure it is right or helps you track down what is wrong.

 

[Bill Simpson]

200 lines in input, 200 lines in output.

Since you appear to be doing something that you might think of as a "moving average" on pairs of points with 200 input points you are only going to get 199 output points.

Or is what you really want for the 200th output the "moving average" of point #200 and point #1? If so then perhaps an appropriately placed If[i<200,stuff,otherstuff] might be of interest.

 

[MartinV]

Ah, calling for a 201th line of array. No wonder it didn't work. I wrote it like this instead:

surf[x_] :=
Do[k = Append[k,
x[[i, 1]]*(x[[i + 1, 2]] - x[[i, 2]]) -
x[[i, 2]]*(x[[i + 1, 1]] - x[[i, 1]])], {i,1,Length[x] - 1},1]; k =
Append[k,
x[[Length[x], 1]]*(x[[1, 2]] - x[[Length[x], 2]]) -
x[[Length[x], 2]]*(x[[1, 1]] - x[[Length[x], 1]])];

I am attempting to create an array "k" that would store all the partial results one after the other. I used Append 199 times (if Length[x] is 200). Then I wrote another line that complets the circle by calculating the 200th and the 1st line and add it at the end of the array "k".

When I evaluate this definition, I get:

Part::partd: Part specification x[[0,1]] is longer than depth of object. >>

Part::partd: Part specification x[[1,2]] is longer than depth of object. >>

Part::partd: Part specification x[[0,2]] is longer than depth of object. >>

General::stop: Further output of Part::partd will be suppressed during this calculation. >>

I guess I'm still doing something wrong.

 

[Bill Simpson]

First, your
Do[blahblahblah, {i,1,Length[x] - 1},1];
I cannot understand what that final ,1 is doing there and I suspect Mathematica cannot either. I assume that is just a typo and you need to remove that ,1. But it excellent that you appear to have posted all the actual code you are using along with the error messages, otherwise this would have been impossible to catch. But you didn't show any error message that would have resulted from that ,1 and I suspect that is again not showing enough to really tell what is going on.

Next, the precedence of ; and function definition I believe is making what you wrote something very different from what you might think. I do not believe both your Do and your final k=blahblahblah; are both part of your surf[x_] function. I believe your surf[x_] function ends at that ; I believe you are going to need () to override the precedence and get these to group the way I think you want. You can, for example, put a Print["oops"]; between those and see if it prints only when you use surf[yourdata] or if it prints the moment you define surf[x_]. If it does the latter that confirms it is not really a part of your surf[x_] definition.

Next, you haven't shown that k needs to be initialized, perhaps to {} or you will get, at least I get, more errors.

If all my guesswork about what I can't see that you have is correct then killing the kernel and restarting the calculation with my own dummy data I see the following:

In[1]:= yy=Table[{j,j+1},{j,1,200}];
surf[x_]:=
(k={};
Do[k=Append[k,x[[i,1]]*(x[[i+1,2]]-x[[i,2]])-x[[i,2]]*(x[[i+1,1]]-x[[i,1]])],{i,1,Length[x]-1}];
k=Append[k,x[[Length[x],1]]*(x[[1,2]]-x[[Length[x],2]])-x[[Length[x],2]]*(x[[1,1]]-x[[Length[x],1]])]
);
surf[yy]

Out[3]= {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,199}

Now my dummy data is completely meaningless, but I include it to show everything that went into and came out of the calculation.

You might also find Reap and Sow interesting tools to master.

yy=Table[{j,j+1},{j,1,200}];
surf[x_]:=
Reap[
Do[Sow[x[[i,1]]*(x[[i+1,2]]-x[[i,2]])-x[[i,2]]*(x[[i+1,1]]-x[[i,1]])],{i,1,Length[x]-1}];
Sow[x[[Length[x],1]]*(x[[1,2]]-x[[Length[x],2]])-x[[Length[x],2]]*(x[[1,1]]-x[[Length[x],1]])]
][[2,1]];
surf[yy]

This will give identical output but without needing your k to accumulate results and the quadratic slowdown associated with Append when you do this lots of times and in this case it would have accidentally avoided your precedence error in your function definition too.

Another technique you might consider when doing things like your "rolling average"

yy = Table[{j, j + 1}, {j, 1, 200}];
surf[x_] :=
Reap[
Do[Sow[x[[i, 1]]*(x[[Mod[i+1,Length[x],1], 2]] - x[[i, 2]]) - x[[i, 2]]*(x[[Mod[i+1,Length[x],1], 1]] - x[[i, 1]])], {i, 1, Length[x]}];
][[2, 1]];
surf[yy]

which you need to study to see all the small changes and test carefully to ensure it always gives exactly the same result.

You could also consider appending the first element of x to the end of x and then doing your "moving average" 200 times. But you will need to watch out when you do this and understand how function arguments really behave in Mathematica, they are not really "variables" that you can then assign new values to inside a function definition, but that a much deeper swimming pool than I want to try to explain at the moment.

All these are opportunities for you to consider how best to use data structures and algorithms to get your result correctly.

 

[MartinV]

I cannot understand what that final ,1 is doing there and I suspect Mathematica cannot either.

It's the unit of the count for i. I looked into Mathematica documentation and I thought I should add it to make sure it goes the way it was suppose to. Not necessary though.

I used the () and I think it's working now. Thanks for the help and your continuous patience. I keep mixing syntax from Matlab because I spent the most time using that.

 

[Bill Simpson]

No.

In[1]:= Do[Print,{i,1,3},1]
From In[1]:= Do::itform : Argument 1 at position 3 does not have the correct form for an iterator.
Out[1]= Do[Print,{i,1,3},1]

That is not going to work, not unless you are not really showing what you actually have, which is a very very bad thing.

But either of these will work.

In[2]:= Do[Print,{i,1,3,1}]
From In[2]:= 1
From In[2]:= 2
From In[2]:= 3

In[3]:= Do[Print,{i,1,3}]
From In[3]:= 1
From In[3]:= 2
From In[3]:= 3

Stepsize is always optional and if not present defaults to 1

 

[MartinV]

OK, this is just unsettling. The code worked fine yesterday. I run it again today and I get this:

ListPlot::lpn: {P,{{-10.045,<<20>>},{<<1>>},<<8>>,<<190>>},<<1>>,{{0.0474999,-0.000911766},<<9>>,<<190>>}} is not a list of numbers or pairs of numbers. >>

Not a list of numbers? Of course it is, I checked. Could Mathematica make mistakes because the code didn't change from yesterday to today.

 

[Dale]

Sometimes when you are testing and changing code you can set some variables and then delete the line where you set them. For instance, here P is is obviously not a number, but perhaps yesterday you had evaluated some statement like:
P=5
So that yesterday P did evaluate to a number, but today it does not.

It is good to quit the evaluation kernel and test your code "fresh" every once in a while.

 

[MartinV]

I found the problem. A typo.

I made a successful animation with additional images appearing and vanishing, just like they should.

All that remains is to use the surf method and see how the surface of those images varies with time.

Thank you for all the help. I could not do it without your assistance.