NDSolve with vectors
Mathematica doesn't have vector variables (yet). That is to say, you can assign a list to a variable, but you cannot use a variable in a function like NDSolve
and let Mathematica work out its dimensions or let the dimensions be undetermined.
If you change your function to this:
gravity = 10;
withDrag[p0_, v0_, drag_] :=
Module[{p},
p[t_] := {p1[t], p2[t], p3[t]};
p[t] /.
NDSolve[
Thread /@ {
p[0] == p0,
p'[0] == v0,
p''[t] == drag*Norm[p'[t]]*p'[t] + {0, 0, -gravity}} // Flatten,
p[t],
{t, 0, 5}
]// First
]
it works. What is does is defining your p
as a vector (list) of functions. Thread
takes care of distributing ==
over the vector components and Flatten
makes a single list of equations from all this.
track[t_] = withDrag[{0, 0, 0}, {0, 10^2, 10}, 0.001];
ParametricPlot3D[track[t], {t, 0, 5}, BoxRatios -> 1]
Note that I reduced the starting value of v0[[2]] to 10^2 because 10^4 yields a 'stiff' system. Also note that I used BoxRatios -> 1
to prevent the box from becoming flat.
While under the hood this method still provides Mathematica with the 9 equations that you already tried manually, it has the advantage that it leaves your vector equations intact.
As of Version 9, you can work with vectors in NDSolve[]
!:
gravity = 10;
withDrag[p0_, v0_, drag_] := Module[{p},
p[t_] := Evaluate@Array[Unique[][t] &, 3];
p[t] /. NDSolve[{
p[0] == p0,
p'[0] == v0,
p''[t] == drag*Norm[p'[t]]*p'[t] + {0, 0, -gravity}},
p[t], {t, 0, 5}] // First]
track[t_] = withDrag[{0, 0, 0}, {0, 10^2, 10}, 0.001];
ParametricPlot3D[track[t], {t, 0, 5}, BoxRatios -> 1]
Having a helper function rhs
, which evaluates only with a numeric vector as argument, for the right-hand side of the force equation lets you use vectors as you want. This way the undesired symbolic precalculation (threading of drag (v.v) Normalize[v]
with {0, 0, gravity}
) is bypassed and the solving continues numerically. See this answer for a bit more detail.
Physically, the drag term should be negative. Also, just as an interesting angle, I added WhenEvent
"equation" that terminates the integration when particle hits the ground.
withDrag[p0_, v0_, drag_] :=
Module[{gravity = 10, rhs},
rhs[v_?(VectorQ[#, NumericQ] &)] :=
-drag (v.v) Normalize[v] - {0, 0, gravity};
NDSolveValue[{
p''[t] == rhs[p'[t]],
p'[0] == v0,
p[0] == p0,
WhenEvent[p[t][[3]] < 0, "StopIntegration"]},
p, {t, 0, \[Infinity]}]]
The solution time depends on the initial values, it can be extracted with suitable parting.
sol = withDrag[{0, 0, 0}, {10, 10, 100}, .1];
ParametricPlot3D[sol[t], {t, 0, sol[[1, 1, 2]]},
BoxRatios -> 1,
ImageSize -> Small]