Damped Orbit

If we introduce polar coordinates and consider the radial and transverse equations of motion, we obtain the equations $\ddot{r} - r\dot{\theta}^2 \; = \; -10r^{-2} - \tfrac{1}{10}\dot{r} \hspace{2cm} 2\dot{r}\dot{\theta} + r\ddot{\theta} \; = \; -\tfrac{1}{10}r\dot{\theta}$ The second of these simplifies to read $\frac{d}{dt}\big(r^2\dot{\theta}e^{\frac{1}{10}t}\big) \; = \; 0$ and hence, given then initial conditions, we deduce that $r^2\dot{\theta} \; = \; \sqrt{10}e^{-\frac{1}{10}t}$ so we see that the angular momentum of the system decays exponentially with time. We can now substitute for $\dot{\theta}$ in the radial equation, obtaining $\ddot{r} - 10r^{-3}e^{-\frac15t} + 10r^{-2} + \tfrac{1}{10}\dot{r} \; = \; 0$ If we attempt the usual substitution $u = r^{-1}$ at this stage, we obtain the equation $\frac{d^2u}{d\theta^2} +u \; = \; e^{\frac15t}$ This is elegant, but not easy to solve, so we revert to solving the differentiaL equation $\ddot{r} - 10r^{-3}e^{-\frac15t} + 10r^{-2} + \tfrac{1}{10}\dot{r} \; = \; 0$ together with the initial conditions $r(0) = 1$ and $\dot{r}(0) = 0$ , numerically. A plot of the first $8$ seconds is below:

There are three times in the first $8$ seconds when the radius reaches the value of $0.2$ , namely $t = \boxed{7.803336}$ , $t = 7.843896$ and $t = 7.972016$ .

This problem was fun to solve. The method using which I obtained the answer can be studied in a commented script of code attached below. Also attached are some plots:

Compare the above to the case where no damping is present:

So basically, the damping results in a reduction of speed which results in a progressive reduction of orbital radius with time. In the absence of damping, we get the classical solution of the two body problem. The script used to generate plots is different from that used to find the result of this problem. In the latter, I ran a simulation until $t=12$ .

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

clear all
clc

% Parameters:
M          = 10;
m          = 1;
G          = 1;

% Initial conditions - Position coordinates and position vector:
x(1)       = 1;
y(1)       = 0;
R(:,1)     = [x(1);y(1);0];

% Initial conditions - Velocity components:
dx(1)      = 0;
dy(1)      = sqrt(10);

% Time, time index and time step initialisation:
t(1)       = 0;
k          = 1;
dt         = 1e-6;


% Running simulation until the distance from the origin is within a specified tolerance:

tol        = 1e-5;

while abs(norm(R(:,k)) - 0.2) >= tol

    Fg         = -(G*M*m)*(R(:,k)/norm(R(:,k))^3);   % Gravity force
    Fd         = -0.1*[dx(k);dy(k);0];               % Damping force

    ddR        = (Fg + Fd)/m;                        % Newton's second law

    % Numerical integration: Semi implicit Euler:

    dx(k+1)    = dx(k) + ddR(1)*dt;
    dy(k+1)    = dy(k) + ddR(2)*dt;

    x(k+1)     = x(k) + dx(k+1)*dt;
    y(k+1)     = y(k) + dy(k+1)*dt;

    R(:,k+1)   = [x(k+1);y(k+1);0];

    t(k+1)     = t(k) + dt;
    k          = k + 1;

end

ANSWER     = t(end);
% ANSWER   = 7.8033

The system of three equations to solve, where r(t) is the distance of the orbiting particle from the origin:

$G = 1$ , $m_1=1$ , $m_2=10$

$m_1 x''(t)=-\frac{G \left(m_1 m_2 x(t)\right)}{\left(x(t)^2+y(t)^2\right)^{3/2}}-0.1 x'(t)$

$m_1 y''(t)=-\frac{G \left(m_1 m_2 y(t)\right)}{\left(x(t)^2+y(t)^2\right)^{3/2}}-0.1 y'(t)$

$r(t)=\sqrt{x(t)^2+y(t)^2}$

Solved it digitally using Wolfram. Here is a plot of the solution for x(t), y(t) and r(t) vs time: