16
$\begingroup$

I'm stumped. I'm trying to write this using vectors, but the 2nd derivative isn't being expanded like I expected it to be. This is a system of equations for a projectile with quadratic drag and gravity (the linear drag is ignored for now). Negative Z is down, X and Y are the horizontal plane. If I write it as 9 equations, one for each coordinate, it works fine, but I'd rather use vectors since it is shorter and (to at least me) more obvious what is going on. Plus since I am new to Mathematica it would be good to learn more/better ways to use it.

gravity = 10; withDrag[p0_, v0_, drag_] := NDSolve[{ p[0] == p0, p'[0] == v0, p''[t] == drag * Norm[p'[t]] * p'[t] + {0,0,-gravity}}, {p}, {t, 0, 5}] withDrag[{0,0,0}, {0,10^4,10}, 0.001]

I get:

NDSolve::ndfdmc: Computed derivatives do not have dimensionality consistent with the initial conditions. >> NDSolve[{ p[0] == {0, 0, 0}, p'[0] == {0, 10000, 10}, p''[t] == { 0.001 Norm[p'[t]] p'[t], 0.001 Norm[p'[t]] p'[t], -10 + 0.001 Norm[p'[t]] p'[t]}}, {p}, {t,0,5}]

I formatted the output to make the error more obvious. Each of elements of the p'' vector has all three elements of p'[t]. Each one should really be p'[t][[dim]] (or something like that).

Any clues as to what I'm doing wrong?

Improve this question
$\endgroup$

4 Answers 4

Reset to default
16
$\begingroup$

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]

Mathematica graphics

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.

Improve this answer
$\endgroup$
6
  • 4
    $\begingroup$ This is also how I would have done it, but it should probably be pointed out that Mathematica does know how to deal with vector functions in some cases. See e.g. this answer. $\endgroup$ Commented Nov 11, 2012 at 22:04
  • $\begingroup$ @jens You're right. I suppose the problem here lies in the assignments with p0 and v0, which aren't explicitly vectors, right? $\endgroup$ Commented Nov 11, 2012 at 22:27
  • $\begingroup$ @jens Vector equations seem to work only if the initial conditions are specified as a scalar constant, not a vector constant. Replace in the doc example the zero in the first example by {0,0,0,0} (which would seem to make more sense) and it fails. $\endgroup$ Commented Nov 11, 2012 at 22:43
  • $\begingroup$ Yes, I guess one could change the function argument from p0_ to {p0x_, p0y_, p0z_} etc., but it seems that even then the second-order differential equation is too hard to recognize as vectorial. So your approach is just the safest, I think. $\endgroup$ Commented Nov 11, 2012 at 22:45
  • $\begingroup$ Nice. I'll have to study the answer some more though... Thanks! $\endgroup$ Commented Nov 12, 2012 at 21:37
16
$\begingroup$

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]

Mathematica graphics

Improve this answer
$\endgroup$
2
$\begingroup$

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]

enter image description here

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ I would advise using explicit variables instead of % (Out) as it might store a different thing on someone else's machine if cells are not evaluated in the exact order of your post. $\endgroup$ Commented Jul 6, 2016 at 7:23
  • $\begingroup$ Can you explain how this works? sol[[1, 1, 2]]? Is this extracting the time it hit ground? $\endgroup$ Commented Mar 22, 2022 at 1:40
2
$\begingroup$

Update

It's a bit surprising that NonThreadable is not needed for this specific case. At least since version 12.3, NDSolve can parse the system in the expected way even if Listable attribute has actually ruined the system:

gravity = 10; withDrag[p0_, v0_, drag_] := NDSolveValue[{p[0] == p0, p'[0] == v0, p''[t] == drag*Norm[p'[t]]*p'[t] + {0, 0, -gravity}}, p, {t, 0, 5}] solp = withDrag[{0, 0, 0}, {0, 100, 10}, 0.001]; ParametricPlot3D[solp[t], {t, 0, 5}, BoxRatios -> 1]

enter image description here

But still, I believe setting a NonThreadable attribute will be a safer choice.

Aha, something cool has been introduced in version 14.1. Now we only need a NonThreadable attribute!:

gravity = 10; SetAttributes[p, NonThreadable] withDrag[p0_, v0_, drag_] := NDSolveValue[{p[0] == p0, p'[0] == v0, p''[t] == drag*Norm[p'[t]]*p'[t] + {0, 0, -gravity}}, p, {t, 0, 5}] solp = withDrag[{0, 0, 0}, {0, 100, 10}, 0.001]; ParametricPlot3D[solp[t], {t, 0, 5}, BoxRatios -> 1]

enter image description here

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.