# InvariantErrorPlot help!

Currently trying to utilise InvariantErrorPlot and failing miserably.

My code is as follows:

Needs["DifferentialEquationsNDSolveUtilities"];
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.8, y1 == 0.6, z1 == 0, x1' == 0,
y1' == 1, z1' == 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]

ExperimentalNumericalFunction::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.}.

ExperimentalNumericalFunction::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.}.

ExperimentalNumericalFunction::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 ExperimentalNumericalFunction::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.

• Is InvariantErrorPlot part of a package? – Mike Honeychurch May 9 '17 at 8:07
• @MikeHoneychurch I've updated the code with the package. It is Needs["DifferentialEquationsNDSolveUtilities"]; – Rumplestillskin May 9 '17 at 9:38

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.8, 0.6, 0},
satpos' == {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'}]
(*  {2, 3, 4, 5, 6, 7}  *)

(* OR *)

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

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


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 =

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

• Well using  Plot[energy /. satorbit, {t, 0, 1000}] seems to work quite well but for some reason the InvariantErrorPlot isn't doing what it's meant to do! hmmmm.. – Rumplestillskin May 9 '17 at 7:22