Damped Oscillation 12-19-2020 (Part 2)

At a general time $t$ , let the coordinates of the particle be $(x,y)$ and the velocity components be $(\dot{x},\dot{y})$ . The equations of motion are:

$m \left[\begin{matrix} \ddot{x} \\ \ddot{y} \end{matrix} \right] = \vec{F}_S + \vec{F}_D + \vec{F}_G$ $\implies m \left[\begin{matrix} \ddot{x} \\ \ddot{y} \end{matrix} \right] = -\frac{K(\sqrt{x^2+y^2}-L_o)}{\sqrt{x^2+y^2} }\left[\begin{matrix} x \\ y\end{matrix} \right] -D\left[\begin{matrix} \dot{x} \\ \dot{y} \end{matrix} \right] +\left[\begin{matrix} 0\\ -mg \end{matrix} \right]$

$x(0)=1 \ ; \ y(0) = 0 \ ; \ \dot{x}(0)=\dot{y}(0)=10$

Solving this numerically for 6.6 seconds yields the required answer which is:

$\mathrm{ANSWER} = \sqrt{x(6.6)^2 + y(6.6)^2}\approx 1.971$

The plot of the particle's trajectory is shown as follows:

One can see that damping makes the motion pretty uninteresting and predictable. The particle eventually settles to its stable equilibrium position of $(0,-3)$ . This is unlike the case of no damping which results is a very nice looking quasi-periodic like trajectory.

I used Newton's laws like Karan Chatrath. The second graph of his is a cool result that is typical of a spring pendulum.

In this scenario, representing the equation of motion in vector form is very useful.

Here are the details of my simulation solution:

Spring force

The spring force acts radially inward towards $(0, 0)$ , so it is useful to develop an expression for the unit radial vector:

$\displaystyle \hat{R} = \braket{\frac{x}{\sqrt{x^2+y^2}}, \frac{y}{\sqrt{x^2+y^2}}}$

So that at any point in the trajectory of the particle we know the direction that the spring force will act in.

As for the final value of the spring force:

$\displaystyle \vec{\bold{F}_s} = -k(\sqrt{x^2+y^2} - l_{0})\hat{R}$

Gravitational force

This one's fairly easy as it's always pointing downwards in our standard basis: $\displaystyle \begin{bmatrix} 0 \\ -mg \\ \end{bmatrix}$

Drag force

The drag force acts in the opposite direction to velocity and has a constant of $0.5$ , so:

$\displaystyle \vec{\bold{F}}_d = -D \begin{bmatrix} \dot{x} \\ \dot{y} \\ \end{bmatrix}$

Final expression

$m \begin{bmatrix} \ddot{x} \\ \ddot{y} \\ \end{bmatrix} = -k(\sqrt{x^2+y^2} - l_{0})\hat{R} + \begin{bmatrix} 0 \\ -mg \\ \end{bmatrix} + -D \begin{bmatrix} \dot{x} \\ \dot{y} \\ \end{bmatrix}$

Now that we are keeping track of $x$ and $y$ and their derivatives, along with useful vectors, we just have to numerically solve, numerically integrating both components and recombining them into a displacement vector.

Here's the code:

 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

time = 0;
deltaT = 10^-6;

%initial conditions
x = 1;
y = 0;

xDot = 10;
yDot = 10;

%constants
m = 1;
g = 10;
k = 5;
l = 1;
D = 0.5;

while time <= 6.6
    %unit radial vector for spring force and displacement
    displacement = [x y];
    unitRadialVector = [x/hypot(x,y) y/hypot(x,y)];

    %forces
    springForce = -k*(hypot(x,y)-l)*unitRadialVector;
    gravitationalForce = [0 -m*g];
    dragForce = -D*[xDot yDot];

    %second law computation
    totalForce = springForce + gravitationalForce + dragForce;

    %numerical integration; second derivatives 
    accelerationVector = totalForce/m;
    xDotDot = accelerationVector(1);
    yDotDot = accelerationVector(2);

    xDot = xDot + xDotDot*deltaT;
    yDot = yDot + yDotDot*deltaT;

    x = x + xDot*deltaT;
    y = y + yDot*deltaT;

    time = time + deltaT;

end

disp(hypot(displacement(1),displacement(2)));