InvariantErrorPlot help!

Question

Currently trying to utilise InvariantErrorPlot and failing miserably.

My code is as follows:

Needs["DifferentialEquations`NDSolveUtilities`"];
T = {t, 0, 1000};
satpos[t_] = {x1[t], y1[t], z1[t]};
satorbit = 
  NDSolve[ Thread /@ {satpos''[t] == -(satpos[t])/Norm[satpos[t] ]^3,
      x1[0] == 0.8, y1[0] == 0.6, z1[0] == 0, x1'[0] == 0, 
      y1'[0] == 1, z1'[0] == 0.5} // Flatten, {x1, y1, z1}, 
   T];
energy = 1/2 Norm[satpos'[t]]^2 - 1/Norm[satpos[t]];
angmom[t] = Norm[Cross[satpos[t], satpos'[t]]];
InvariantErrorPlot[energy, {x1[t], y1[t], z1[t]}, T, satorbit, 
 PlotStyle -> Green]

Im not getting (any) desired output. Any suggestions on a quick fix?

The output is given below

InvariantErrorPlot[energy, satpos[t], t, satorbit]

Experimental`NumericalFunction::nlnum: The function value {-1.+1/2 (Abs[(x1^\[Prime])[0.]]^2+Abs[(y1^\[Prime])[0.]]^2+Abs[(z1^\[Prime])[0.]]^2)} is not a list of numbers with dimensions {1} at {t,x1[t],y1[t],z1[t]} = {0.,0.8,0.6,0.}.

Experimental`NumericalFunction::nlnum: The function value {-1.+1/2 (Abs[(x1^\[Prime])[0.]]^2+Abs[(y1^\[Prime])[0.]]^2+Abs[(z1^\[Prime])[0.]]^2)} is not a list of numbers with dimensions {1} at {t,x1[t],y1[t],z1[t]} = {0.,0.8,0.6,0.}.

Experimental`NumericalFunction::nlnum: The function value {-0.978827+1/2 (Abs[(x1^\[Prime])[0.0361774]]^2+Abs[(y1^\[Prime])[0.0361774]]^2+Abs[(z1^\[Prime])[0.0361774]]^2)} is not a list of numbers with dimensions {1} at {t,x1[t],y1[t],z1[t]} = {0.0361774,0.799488,0.635785,0.0180849}.

General::stop: Further output of Experimental`NumericalFunction::nlnum will be suppressed during this calculation.

Transpose::nmtx: The first two levels of {0,<<49>>,<<1180>>} cannot be transposed.

Following this, there is a graph which looks quite correct but I'm not following the output messages. Any ideas?

Edit simply using Plot[energy/.satorbit, {t,0,100}] works well.

Answer 1

Take Akku14's suggestion a bit further: In a second-order ODE, one sometimes imagines the first derivatives are variables. Apparently you have to tell InvariantErrorPlot and NDSolve explicitly all variables:

Clear[t, x1, y1, z1];
satorbit = NDSolve[        (* refactored *)
   {satpos''[t] == -(satpos[t])/Norm[satpos[t]]^3,
    satpos[0] == {0.8, 0.6, 0},
    satpos'[0] == {0, 1, 0.5}},
   {x1, y1, z1, x1', y1', z1'}, {t, 0, 1000}];

InvariantErrorPlot[energy, {x1, y1, z1, x1', y1', z1'}, t, satorbit, PlotStyle -> Green]

Irritatingly InvariantErrorPlot assigns values to the variables, which should be considered a bug:

Through[{x1, y1, z1, x1', y1', z1'}[1]]
(*  {2, 3, 4, 5, 6, 7}  *)

(* OR *)

DownValues /@ {x1, y1, z1}
(*  {{HoldPattern[x1[1]] :> 2}, {HoldPattern[y1[1]] :> 3}, {HoldPattern[z1[1]] :> 4}}  *)

SubValues@Derivative
(*
  {HoldPattern[Derivative[1][x1][1]] :> 5, 
   HoldPattern[Derivative[1][y1][1]] :> 6, 
   HoldPattern[Derivative[1][z1][1]] :> 7}
*)

Answer 2

Since in your "energy" the derivatives x1'[t], y1'[t], z1'[t] are used, you have to let NDSolve give pure functions x1,y1,z1 as answer.

satorbit = 
   NDSolve[Thread /@ {satpos''[t] == -(satpos[t])/Norm[satpos[t]]^3, 
  x1[0] == 0.8, y1[0] == 0.6, z1[0] == 0, x1'[0] == 0, 
  y1'[0] == 1, z1'[0] == 0.5} // Flatten, {x1, y1, z1}, T, 
    MaxSteps -> 20000];

(Appendix: I could'nt test with "InvariantErrorPlot", but only with Plot, since I used MMA Version 8.0)

 Plot[energy /. satorbit, {t, 0, 1000}]