# How to solve $\dot{\mathbf{v}} = - \frac{GM}{r^3}\mathbf{r}$ using differential eq-n solver in MATLAB

I need to solve second order system with boundary conditions with help of differential equation solver. The system is $$\dot{\mathbf{v}} = - \frac{GM}{r^3}\mathbf{r}.$$ I tried to use MATLAB dsolve command like this:

syms x(t) y(t)
eqns = [diff(x,t,2) == -1*x/sqrt(x*x+y*y)^3, diff(y,t,2) == -1*y/sqrt(x*x+y*y)^3];
Dx = diff(x,t);
Dy = diff(y,t);
cond = [x(0) == 10, y(0) == 0, Dx(0) == 0, Dy(0) == 1];
dsolve(eqns,cond)
Warning: Explicit solution could not be found.


EDIT: I think i need to use ode45 instead, something like this:

f = @(t,y)[y(3); y(4); -y(1)/sqrt(y(1)*y(1)+y(2)*y(2))^3; -y(2)/sqrt(y(1)*y(1)+y(2)*y(2))^3];
[ts,ys]=ode45(f, [0, 100], [10; 0; 0; 0.1]);
plot(ys(:,1), ys(:,2));


And it kinda works:

Very good tutorial which helped me a lot: Using Matlab for Higher Order ODEs and Systems of ODEs

• Do you want to calculate the two body problem? Commented Mar 6, 2017 at 16:19
• @MrYouMath yes, you are right. If you have anything interesting to add, I will be glad to upvote and accept your answer. Thanks!
– Yola
Commented Mar 6, 2017 at 16:35
• To test your code: the initial condition you have used here corresponds to the analytical solution $x = 10\cos\left(\frac{t}{100}\right)$ and $y = 10\sin\left(\frac{t}{100}\right)$. Commented Mar 6, 2017 at 17:55
• @Winther Your solution looks like a circle for me, not sure it is right.
– Yola
Commented Mar 7, 2017 at 8:46
• Yes it corresponds to slightly different initial conditions. In general if you fix $R$ and take $\omega = \frac{1}{R^{3/2}}$ and take the initial conditions to be $[R,0,0,R\omega]$ then the analytical solution is $x(t) = R\cos(\omega t)$, $y(t) = R\sin(\omega t)$. You need to evolve it over $t = 0$ to $t = \frac{2\pi}{\omega}$ to get one full orbit. btw you should increase the accuracy when solving here. Try something like options = odeset('RelTol',1e-8,'AbsTol',1e-8); and append options to the call as ode45( .... , options). Commented Mar 7, 2017 at 17:27

You can use

function ydot = orbit(~, y)
%ORBIT: differential equation for 2D orbit
%  Returns the vector of the derivatives of y at time t.
ydot = [y(3:4); -1/norm(y(1:2)) * y(1:2)];
end


When you call ode45 (with @orbit) you may want to set the MaxStep option (with odeset) to get better accuracy. Of course, you can also write an anonymous function directly.

EDIT: I just saw your edit. The non-closed orbit is probably the result of not controlling the maximum integration step.

Another way is building a matrix $3 \times N$, and a digital differentiation filter

$$\frac{1}{\Delta t}\begin{bmatrix}1&-1\end{bmatrix}$$

Denote the matrix representation $\bf D$ of it's convolution on a vectorization of said image $\bf v$. To be strict, we should also use some interpolation filter to make sure we measure differential and function on the same points in time:

$$\frac{1}{2}\begin{bmatrix}1&1\end{bmatrix}$$

Let us call the matrix representation of this interpolation $\bf L$. Now let us solve a sequence of norm problems. Norm 2 (resulting in least-squares) should probably be good enough for just a simple test. $${\bf v_0} = \min_{\bf v}\left\{\left\|\frac{\bf Lv}{d^3}-\bf Av\right\| + {\bf C_i}\|\bf v - v_i\|\right\}$$ Where $d^3 = |{\bf v}|^3$ updated after each new solved one. $\bf C_i$ is a matrix encoding the certainty of boundary (initial value) condition and $\bf v_i$ the values of said initial/boundary values.

This equation system will have a set of equations "grabbing ahold" of neighbouring time points implicitly giving a very powerful processing.

I'll get back if I get a chance to try it.