3
$\begingroup$

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.

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

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]

Mathematica graphics

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}
*)
$\endgroup$
2
$\begingroup$

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}]
$\endgroup$
  • $\begingroup$ 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.. $\endgroup$ – Rumplestillskin May 9 '17 at 7:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.