# Interpolating initial conditions of orbits

I have a data file containing the $(x,y,p_x,p_y)$ initial conditions of two-dimensional orbits.

Let's plot them on the $(x,y)$ plane.

data0 = Import["nhim2.dat", "Table"];
d0 = data0[[All, {1, 2}]];
S0 = ListPlot[d0, PlotStyle -> {Blue, PointSize[0.003]}]


Then we define some random $(x_0,y_0)$ initial conditions on the same plane

data = {{0.9896, -0.29035}, {0.9214, 0.20255}, {0.96015, 0.39785},
{0.96325, 0.325}, {1.0764, 0.0367}, {1.00045, 0.2227},
{0.8253, -0.3167}, {1.0361, 0.36065}, {0.83925, -0.02685},
{1.02525, -0.34305}, {1.0113, 0.10645}, {1.0082, -0.1431},
{0.8284, 0.22115}, {0.90435, 0.13435}, {0.90125, -0.1369},
{1.0206, 0.2475}, {1.1198, -0.02685}, {0.816, 0.2909},
{0.8966, 0.3839}, {0.8129, 0.1049}};

S1 = ListPlot[data, PlotStyle -> {Red, PointSize[0.015]}];
Show[{S0, S1}]


My question is the following: How can I use the initial data0 as a base in order to interpolate the $p_x$ and $p_y$ of the 20 random $(x_0,y_0)$ initial conditions? The aim is to create a new list, data2, containing in four columns the $(x,y,p_x,p_y)$ of the 20 points.

I use version 9 of Mathematica in Win XP SP3.

UPDATE

The answer proposed by @JasonB seems to work. However for this data file there are some issues.

Now the initial conditions on the $(x,y)$ look like

data0 = Import["nhim33.dat","Table"];
d1 = GatherBy[data0, #[[;; 2]] &] // GatherBy[#, Length] & // Last //
Flatten[#, 1] &;
Length[d1]


I end up with only 48 points from the initial 27407.

Then the

pyfunc = Interpolation[{{#1, #2}, #4} & @@@ d1, InterpolationOrder -> 1];
pxfunc = Interpolation[{{#1, #2}, #3} & @@@ d1, InterpolationOrder -> 1];


cannot be evaluated without errors.

Any ideas?

• Is it possible to easily split up data0 into sublists for the individual orbits? – Jason B. Jun 23 '16 at 16:01
• @JasonB Unfortunately no. data0 is just a bunch of initial conditions. In this case, they form nice ellipses but in other case they are randomly scattered points. If we use Interpolation[data0] would it work? – Vaggelis_Z Jun 23 '16 at 16:03
• Then I'm confused how this differs from the previous question (I had thought the ellipses were orbits and you wanted to interpolate orbits). Just from reading, I would recommend the same strategy as in that question - make an interpolating function on the unstructured grid given by d0. Also, your data0 doesn't have the px or py values. – Jason B. Jun 23 '16 at 16:08
• @JasonB I just uploaded the correct file. The interpolation should take into account the px and py of data0 and then apply it to the random 20 initial conditions. I think that this is different form the previous case. Anyway, provide an answer so as to see whether we mean the say the same thing! – Vaggelis_Z Jun 23 '16 at 16:16
• @JasonB When I try pyfunc = Interpolation[{{#1, #2}, #4} & @@@ data0, InterpolationOrder -> 1] I get the following error message: "Interpolation::indp: There are duplicated abscissa points in ..." – Vaggelis_Z Jun 23 '16 at 16:22

A couple of points:

1. You have some duplicate coordinates in your data, these need to be removed before you can interpolate
2. A couple of the points in your list are outside the range you've specified. Therefore you need to extrapolate, not interpolate, and the results at these points can't be trusted.

data0 = DeleteDuplicatesBy[data0, #[[;; 2]] &];
pyfunc = Interpolation[{{#1, #2}, #4} & @@@ data0,
InterpolationOrder -> 1];
pxfunc = Interpolation[{{#1, #2}, #3} & @@@ data0,
InterpolationOrder -> 1];

data2 = {#1, #2, pxfunc[#1, #2],
pyfunc[#1, #2]} & @@@ {{0.9896, -0.29035}, {0.9214,
0.20255}, {0.96015, 0.39785}, {0.96325, 0.325}, {1.0764,
0.0367}, {1.00045, 0.2227}, {0.8253, -0.3167}, {1.0361,
0.36065}, {0.83925, -0.02685}, {1.02525, -0.34305}, {1.0113,
0.10645}, {1.0082, -0.1431}, {0.8284, 0.22115}, {0.90435,
0.13435}, {0.90125, -0.1369}, {1.0206,
0.2475}, {1.1198, -0.02685}, {0.816, 0.2909}, {0.8966,
0.3839}, {0.8129, 0.1049}}
(* {{0.9896, -0.29035, 0.0551768, 0.80115}, {0.9214,
0.20255, -0.0651558, 1.11407}, {0.96015, 0.39785, -0.0828736,
0.83324}, {0.96325, 0.325, -0.0722046, 0.878227}, {1.0764,
0.0367, -0.00358338, 0.546399}, {1.00045, 0.2227, -0.0414146,
0.797068}, {0.8253, -0.3167, 0.152263, 1.40477}, {1.0361,
0.36065, -0.0430032, 0.554697}, {0.83925, -0.02685, 0.0143835,
1.49376}, {1.02525, -0.34305, 0.045252, 0.63029}, {1.0113,
0.10645, -0.0194832, 0.796962}, {1.0082, -0.1431, 0.0264175,
0.799015}, {0.8284, 0.22115, -0.116249, 1.46361}, {0.90435,
0.13435, -0.0493386, 1.21041}, {0.90125, -0.1369, 0.0511158,
1.22164}, {1.0206, 0.2475, -0.0380753, 0.705899}, {1.1198, -0.02685,
0.00127011, 0.370953}, {0.816, 0.2909, -0.147451,
1.44989}, {0.8966, 0.3839, -0.121433, 1.07842}, {0.8129,
0.1049, -0.0633001, 1.58377}} *)


For version 9, here is another way of getting the duplicate-abscissa free version:

data0 = GatherBy[data0, #[[;; 2]] &] // GatherBy[#, Length] & //
Last // Flatten[#, 1] &;

• Unfortunately, DeleteDuplicatesBy is not recognized in v9 – Vaggelis_Z Jun 23 '16 at 16:26
• Try using the much slower two-argument version of DeleteDuplicates - DeleteDuplicates[data0, #1[[;; 2]] == #2[[;; 2]] &] – Jason B. Jun 23 '16 at 16:28
• @Vaggelis_Z - see the edit for a faster version. – Jason B. Jun 23 '16 at 16:34
• When I evaluate data2 = {#1, #2, pxfunc[#1, #2], pyfunc[#1, #2]} & @@@ data I get the error message "InterpolatingFunction::dmval: "Input value {0.9896,-0.29035} lies outside the range of data in the interpolating function. Extrapolation will be used" and then all py form the new list is 0. – Vaggelis_Z Jun 23 '16 at 16:40
• I get that error, but I do get answers as well. I'm using a different version 9 so maybe that's the trouble. Why don't you just remove the points that are outside the original region? – Jason B. Jun 23 '16 at 16:48

This is not a full answer, but it seems to me that you may have a use for code that allows you to group the points by the orbit they belong to, as Jason also mentioned. It seems to me that, with this in hand, the problem may be reduced to the one @JasonB solved in your previous question.

Here's a start in that direction:

paths = FindClusters[data0, 15, Method -> "Agglomerate"];
ListPlot@paths
ListPlot /@ paths