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

Question

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

Answer 1

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.

Answer 2

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.