# Using Grad and NDSolve with vector variables

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]
}
]


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?

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}]]