Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a system of differential equations (referred to as "s") and use NDSolve to obtain the solution. I substitute the interpolated functions for the original functions in order to represent the trajectories of two objects in 3 dimensions. Here is the system that I am working with:

(I added this as an edit after I got some feedback that some more information about what I was attempting to animate would be useful )

sysstars = {((sx[t]^2 + sy[t]^2 + sz[t]^2)^(3/2))*(1 + 
      2*(mstar1/mstar2) + (mstar1/mstar2)^2)*sx''[t] == 
   g*mstar2*(-sx[t]), ((sx[t]^2 + sy[t]^2 + sz[t]^2)^(3/2))*(1 + 
      2*(mstar1/mstar2) + (mstar1/mstar2)^2)*sy''[t] == 
   g*mstar2*(-sy[t]), ((sx[t]^2 + sy[t]^2 + sz[t]^2)^(3/2))*(1 + 
      2*(mstar1/mstar2) + (mstar1/mstar2)^2)*sz''[t] == 
   g*mstar2*(-sz[t]), sx[0] == tpostar1[[1]], sy[0] == tpostar1[[2]], 
  sz[0] == tpostar1[[3]], sx'[0] == v0com1[[1]], 
  sy'[0] == v0com1[[2]], 
  sz'[0] == 
   v0com1[[3]], ((qx[t]^2 + qy[t]^2 + qz[t]^2)^(3/2))*(1 + 
      2*(mstar2/mstar1) + (mstar2/mstar1)^2)*qx''[t] == 
   g*mstar1*(-qx[t]), ((qx[t]^2 + qy[t]^2 + qz[t]^2)^(3/2))*(1 + 
      2*(mstar2/mstar1) + (mstar2/mstar1)^2)*qy''[t] == 
   g*mstar2*(-qy[t]), ((qx[t]^2 + qy[t]^2 + qz[t]^2)^(3/2))*(1 + 
      2*(mstar2/mstar1) + (mstar2/mstar1)^2)*qz''[t] == 
   g*mstar2*(-qz[t]), qx[0] == tpostar2[[1]], qy[0] == tpostar2[[2]], 
  qz[0] == tpostar2[[3]], qx'[0] == v0com2[[1]], 
  qy'[0] == v0com2[[2]], qz'[0] == v0com2[[3]]}

s = NDSolve[
  sysstars, {sx[t], sy[t], sz[t], qx[t], qy[t], qz[t]}, {t, 0, 
   8*10^7}, MaxSteps -> 10^8]  

And here is how I am plotting them:

 Evaluate[{sx[t], sy[t], sz[t], qx[t], qy[t], qz[t]} /. s, {t, 0, 
   10^7}], PlotRange -> All]
 Evaluate[{qx[t] + comx[t], qy[t] + comy[t], qz[t] + comz[t], 
    sx[t] + comx[t], sy[t] + comy[t], sz[t] + comz[t]} /. s, {t, 0, 
   3*10^7}], PlotRange -> All,]

(The second Parametric Plot above is where the center of mass vector is added back to show the system from outside the center of mass inertial frame. So the first one plots two objects on initial trajectories that then go into orbit around the center of mass, with the center of mass treated as static)

The results of this are good - static lines in 3D, with the sensible paths. Now, I would like to animate this so that, say two spheres, move as a function of t. I have tried a few on-line suggestions using Manipulate[Show[..., but without success so far. I believe that the difficulty may relate to the use of Evaluate inside ParametricPlot3D - in conjunction with the animation techniques. It does not seem possible to get rid of Evaluate with ParamatricPlot3D, and still get ParametricPlot3D to work with the interpolated functions which get substituted to generate the representation. There must be a way to do this in Mathematica - the program is way to good not to have a good way to do such an obvious thing with solutions, but I am stumped.

Suggestions greatly appreciated (after too many hours) !


share|improve this question
Welcome to Mathematica.SE! In case you didn't know, you can format your code better by putting four spaces at the front of every code block (you can use the curly brace button above the editing area), or wrapping short code snippets in backticks ``. This will make your post easier to read. For a more extensive introduction to the mark-up system, read this. – rcollyer Aug 5 '12 at 18:52
It's hard to tell without knowing exactly what the issue is - but this could be a duplicate of Animate ParametricPlot3D for two different parametric equations – Jens Aug 5 '12 at 20:04
you can format your code as rcollyer pointed out. I had done this but you reverted it. also, without defining tpostar1 this code doesn't run. finally, it is probably more useful to put a minimal example of what you'd like to do rather than all the actual code you have. – acl Aug 5 '12 at 20:39
Off topic: It might come out a better result if adopting some symplectic methods or invariant preserved methods for this kind of ODE. – Silvia Aug 6 '12 at 13:16


Please note that your eq. could be written in a compact form:

eq[s_List, m1_, m2_, tPos_List, v0Com_List, g_] :=
  With[{st = Through[s[t]]},
   ((Norm[st]^3 (1 + 2 m1/m2 + (m1/m2)^2) #''[t] == -g m2 #[t]) & /@ s)
    Thread[Through[s[0]] == tPos]
    (Thread[D[st, t] == v0Com] /. t -> 0)

And solve them as

s = NDSolve[Join[
    eq[{sx, sy, sz}, 1, 1.7, {1, -1, 1}, {1, 1, 1}, 10],
    eq[{qx, qy, qz}, 1.7, 1, -{1, -1, 1}, -{1, 1, 1}, 10]],
   {sx, sy, sz, qx, qy, qz}, {t, 0, 8}, MaxSteps -> 10^3];

For Example

s = Table[k[j] -> Interpolation[Table[{i, Sin[i j]}, {i, 6}]], {j, 1, 6}]; 
 ParametricPlot3D[{{k[4][t], k[5][t], k[6][t]}, {k[1][t], k[2][t], k[3][t]}} /. s, {t, 1, a}, 
  PlotStyle -> {Red, Green}, 
  Evaluated -> True, 
  PlotRange -> {{-1.5, 1.5}, {-1.5, 1.5}, {-1.5, 1.5}}], {a, 1.01, 6}]

enter image description here


Solving your equations:

(*Some constants*)
g = 10;
mstar1 = 1;
mstar2 = 1.5;
tpostar1 = {1, -1, 1};
tpostar2 = -{1, -1, 1};
v0com2 = {1, 1, 1};
v0com1 = {-1, -1, -1};

Your def for sysstars = ... comes here ....

s = NDSolve[sysstars, {sx[t], sy[t], sz[t], qx[t], qy[t], qz[t]},{t, 0, 8},MaxSteps -> 10^3]

 ParametricPlot3D[{{sx[t], sy[t], sz[t]}, {qx[t], qy[t], qz[t]}} /. s[[1]], {t, 1, a},
  PlotStyle -> {Red, Green},
  PlotRange -> 5 {{-1, 1}, {-1, 1}, {-1, 1}}, 
  Evaluated -> True],
 {a, 0, 10}]

enter image description here


Adding my orbiting balls:

tt = Manipulate[Show@{
    ParametricPlot3D[{{sx[t], sy[t], sz[t]}, {qx[t], qy[t], qz[t]}} /. s[[1]], {t, 1, a},
     PlotStyle -> {Red, Green},
     PlotRange -> 8.5 {{-1, 1}, {-1, 1}, {-1, 1}}, Evaluated -> True],
    Graphics3D[{Sphere[{sx[t], sy[t], sz[t]}], 
                Sphere[{qx[t], qy[t], qz[t]}]} /. s[[1]]] /. t -> a},
  {a, 0, 10}]

enter image description here

share|improve this answer
belisarius - thanks a bunch, for the help – Delphi Aug 5 '12 at 21:30
Really like the orbiting balls...That is really nice touch. I am going to see if I can define some that are similar and put in.... – Delphi Aug 5 '12 at 23:33

Without knowing anything about your functions, I can see that the code you showed isn't quite correct: The Evaluate should wrap only the first plot argument and not the second one containing the plot range. Also, Evaluate can't be followed by { as you typed - it has to be [. Finally, the grouping of the vector components has to be in triplets as shown in the answer by belisarius.

Moreover, when making any kind of trajectory animation you have to fix the PlotRange and AspectRatio of your plot frames to a constant value, otherwise there will be no stationary reference frame to which the movement of your objects can be compared. This means you can't use PlotRange -> All here - except in the first plot over the full time range, to get an idea of the final plot range you're going to need.

After that, to make the animation, you should replace this by a fixed PlotRange as belisarius also shows.

share|improve this answer

Your Answer


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.