1
$\begingroup$

A simple orbital model might be developed like this:

day = 86400;
G=6.67408*10^-11;
ER=6367440;
earthmass=5.9722*10^24;

pot[{x_, y_, z_}] := -G earthmass / Sqrt[x^2 + y^2 + z^2];
grad = -Grad[pot[{x[t], y[t], z[t]}], {x[t], y[t], z[t]}];
odesys = { {x[0], y[0], z[0]} == {1.1 ER, 0, 0}, {x'[0], y'[0], 
     z'[0]} == 
    RotationTransform[ 
      98 \[Degree], {1, 0, 0}]@(7543.7 {0, 1, 0}), {x''[t], y''[t], 
     z''[t]} == grad} ;
soln = NDSolve[odesys, {x, y, z}, {t, 0, 5 day} ];

And that works just fine:

Show[ { 
  Graphics3D[ {Opacity[0.5], Sphere[{0, 0, 0}, ER]}],
  ParametricPlot3D[{x[t], y[t], z[t]} /. soln, {t, 0, 5 day}, 
   PlotRange -> Table[{-2 ER, 2 ER}, 3], AspectRatio -> Automatic]
  }
 ]

gradient orbital

But what I would like to do is just use the variables as vectors, something like:

odesys = {x[0] == {1.1 ER, 0, 0}, x'[0] == {0, 7534, 0}, 
  x''[t] == -earthmass Norm[x[t]]^-2 G Normalize[x[t]]}  

which works

NDSolve[odesys, x, {t, 0, 5 day}]

But I'm not sure how to go about getting Grad to take a variable and make the assumption that it is a 3-element vector. How should I go about this?

$\endgroup$
1
$\begingroup$

Not a perfect answer, but requires less manual analysis compared to yours. Notice that Grad can be replaced with D:

Grad[f[x, y, z], {x, y, z}] == D[f[x, y, z], {{x, y, z}}]
(* True *)

and D works partly correctly when we use an implicit array as 2nd argument:

D[vector.vector, {vector}] /. vector -> {x, y, z}
(* 1.{x, y, z} + {x, y, z}.1 *)

To amend this, we introduce a replace rule:

rule = Dot[a___, 1, b___] :> Dot[a, IdentityMatrix[3], b];

odesys = {x[0] == {1.1 ER, 0, 0}, x'[0] == {0, 7534, 0}, 
   x''[t] == D[G earthmass/Sqrt[x[t].x[t]], {x[t]}]} /. rule

solnvector = NDSolve[odesys, x, {t, 0, 5 day}][[1]]

With[{x = x /. solnvector}, 
 Animate[Show[{Graphics3D[{Opacity[0.5], Sphere[{0, 0, 0}, ER], Red, PointSize@Large, 
        Point@x[t]}], #}], {t, 0, 5 day}] &@ParametricPlot3D[x[t], {t, 0, 5 day}]]

enter image description here

$\endgroup$

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.