I have the following system of ODEs used the propagate the trajectory of a particle in free space about a gravitating body:

$$x'(t) = v_{x}(t)$$ $$y'(t) = v_{y}(t)$$ $$v_{x}'(t) = \frac{-G M x(t)}{(x(t)^{2}+y(t)^{2})^{3/2}}$$ $$v_{y}'(t) = \frac{-G M y(t)}{(x(t)^{2}+y(t)^{2})^{3/2}}$$

Where $t$ is time, $x$ and $y$ are the Cartesian coordinate points of the particle, $v_{x}$ and $v_{y}$ are the $x$- and $y$-velocity components of the particle, $G$ is the universal gravitational constant and $M$ is the mass of the gravitating body (for example, Earth). A 3D example of a particle orbiting a larger mass body is shown below:

Example of a particle orbiting a larger body

If I solve this system of ODEs by integrating the equations using, for example, MATLAB's ode45 solver, I get values for $x$, $y$, $v_{x}$ and $v_{y}$ with respect to time $t$. However, I would also like to be able to find values for the radial velocity, $v_{r}$, of the particle with respect to time $t$ (that is, the speed at which the particle moves towards or away from the gravitating body along a radial line), without having to solve the system of ODEs using a polar coordinate formulation.

enter image description here

As such, is there a way to formulate the radial acceleration, $v_{r}'= dv_{r}/dt$ of the particle, using the available Cartesian states $(x(t),y(t),v_{x}(t),v_{y}(t))$, as a 5th equation to be added to the original system of 4 equations, so that when the system is integrated, I receive output for the radial velocity of the particle as well?

So, for example, the new system of equations becomes:

$$x'(t) = v_{x}(t)$$ $$y'(t) = v_{y}(t)$$ $$v_{x}'(t) = \frac{-G M x(t)}{(x(t)^{2}+y(t)^{2})^{3/2}}$$ $$v_{y}'(t) = \frac{-G M y(t)}{(x(t)^{2}+y(t)^{2})^{3/2}}$$ $$v_{r}'(t) = ?$$

Any help would be much appreciated

EDIT: The following formulation given by nnovich-OK is possible in Mathematica using:

$r'(t) = \frac{d}{dt}\sqrt{x(t)^{2}+y(t)^{2}}$

which is shown in the image below for a snippet of Mathematica code that numerically integrates the system of 5 ODEs:

enter image description here

However, I have not been able to get such a formulation to work in MATLAB, and as such I was hoping for a more "traditional" approach without having to take the derivative of the radius with respect to time.


$v_r$ is tied to $x(t)$ and $y(t)$, which you already found by solving system of ODE: $$ r(t) = \sqrt{x(t)^2 + y(t)^2} $$ $$ v_r (t) = r'(t) = (\sqrt{x(t)^2 + y(t)^2})' $$

| cite | improve this answer | |
  • $\begingroup$ Thanks nnovich-OK. The form $\frac{d}{dt}\sqrt{x(t)^{2}+y(t)^{2}}$ definitely works in Mathematica, however, I have been unable to get it to work in MATLAB (I'll update my original post to give an example) $\endgroup$ – InquisitiveInquirer Apr 12 '17 at 18:36
  • $\begingroup$ Sadly I'm not familiar with MATLAB. If trouble is complexity of the formula, you can try manually simplified variant: $v_r (t) = \frac{x(t) v_x(t) + y(t) v_y(t)}{\sqrt{x(t)^2+y(t)^2}}$ $\endgroup$ – nnovich-OK Apr 12 '17 at 18:45
  • $\begingroup$ The above equation you just gave for $v_{r}(t)$ matches my data well, thanks very much nnovich-OK, problem solved. $\endgroup$ – InquisitiveInquirer Apr 12 '17 at 18:51

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.