I'm trying to use Mathematica to solve a problem of orbits, specifically solving for the trajectory of a spacecraft moving near the Earth and Moon. I assume that the Earth is fixed at the center of the system and the Moon revolves around it, then model the motion of the spacecraft using a vector {x,y,z}. Unfortunately, NDSolve throws an error when I try to account for the gravity caused by the Moon. If I only include the terms for the Earth's gravity, it works perfectly, but it fails otherwise. Below is the code I used:
M[t_] = {rm Cos[((2 Pi)/pm) t], rm Sin[((2 Pi)/pm) t], 0.}
orbit = NDSolve[
{l''[t] == (G me (-l[t]))/(Norm[
l[t]]^3) + (G mm (M[t] - l[t]))/(Norm[M[t] - l[t]]^3),
l[0] == p0,
l'[0] == v0}
, l[t], {t, 0, tmax}, MaxSteps -> 1000000,
Method -> "StiffnessSwitching"];
All of the constants here (rm, pm, G, me, mm, tmax) were defined before running this cell. I'm new to Mathematica and don't quite understand what's going wrong. This is the relevant error:
NDSolve::ndfdmc: Computed derivatives do not have dimensionality consistent with the initial conditions.
M[t_?NumericQ] := {rm Cos[((2 Pi)/pm) t], rm Sin[((2 Pi)/pm) t], 0.}. See https://mathematica.stackexchange.com/questions/18393/what-are-the-most-common-pitfalls-awaiting-new-users/26037#26037 – Michael E2 Jun 07 '20 at 20:21versus what happens if you substitute all at oncea+b /. {a -> {10,20,30}, b -> {1,2,3}}. That's like what's happening in your code, but simplified. It's because{10,20,30}+bis treated as an array operation, "add the scalarbto each element." It doesn't wait to find out ifb` is going to actually be a scalar or a vector. Besides the successive operation is useful sometimes. All this flexibility does create pitfalls for the newbie, though. – Michael E2 Jun 07 '20 at 23:15