#   Numerical investigation of an oscillator
We consider a particle of mass $m = 1$ moving on the horizontal $x$ axis in a potential $U=\omega_{0}^2 x^{2}/2$,
subject to a friction force $F_{friction}=-2\gamma \dot{x}$ and a driving force $F_{driving}=\alpha\cos(t+\phi)$. The equation of motion is therefore of the general form

\begin{equation} \ddot{x}+2\gamma \dot{x}+ \omega_0^2x-\alpha\cos(t+\phi)=0  \qquad(1)\end{equation}

This is a harmonic oscillator with damping (friction) and driving. We want to study the solutions of this equation in terms of the free parameters of the problem. These are the initial position $x_{0}$, the initial velocity $\dot{x}_{0}$, the frequency $\omega_{0}$, and the parameters $\gamma$, $\alpha$, and $\phi$. Once we have the
solution we want to produce plots of the position, the velocity, the potential and kinetic energy,and the phase portrait (velocity in terms of position).
### - First step:
By denoting the velocity $v$ (i.e. $v = \dot{x}$), rewrite (1) as a system of first order ODEs in $x$ and $v$.



$\frac{dx}{dt} = v$ 

$\frac{dv}{dt}$ = -2$\gamma v$ - $w^{2}_{0}x$ + $\alpha\cos(t+\phi)$

### - Second step:  
Let $\{t_k \}$ be a partition of $[0,60]$ such that $0=t_0<t_1<\cdots<t_{600}=60$ and $H$ be the constant length of the $k$-th subinterval (i.e. $H = t_k -t_{k-1}$).

Using <b>odeint</b>, write a Python function named <b>Oscillator_func</b> which takes $x_{0}$, $\dot{x}_{0}$, $\omega_{0}$, $\gamma$, $\alpha$ and  $\phi$ as inputs and return the arrays $x$ and $v$ (numerical solution of the system of ODEs obtained above) such that for a given index $k$ we have $x[k]= x(t_k)$ and $v[k] = v(t_k)$




```python
from scipy.integrate import odeint
from pylab import*
def F (s,t):
    [x,v] =s
    dxdt = v
    dvdt = -2*gamma*v -(w_0)**2*x + alpha*(np.cos(t+phi))
    return [dxdt,dvdt]
def Oscillator_func(x_0,x_dot_0,gamma,w_0,alpha,phi):
    s0 = [x_0,x_dot_0]
    N = (int((b-a)/H))+1
    t=np.linspace(0,60,N)
    sol = odeint(F,s0,t)
    x = sol[:,0]
    v = sol[:,1]
    return x,v,t

```

### - Third step: 
You can then play with the numerical values ($x_{0}$, $\dot{x}_{0}$, $\omega_{0}$, $\gamma$, $\alpha$, $\phi$) to investigate different situations. From now on, we assume that

 \begin{equation} x_{0}=5,\qquad \dot{x}_{0}=0,\qquad \omega_{0}=0.8,\qquad \phi=10.\end{equation} 
 
 
write a code which produces the following three graphs:

 - First graph:
position $x(t)$ and velocity $\dot{x}(t)$, both as a function of time.

 - Second graph: 
 potential energy $U(t)$ and kinetic energy $K(t)$, both as a function of time.
 
 - Third graph: 
 a parametric plot representing the position $x(t)$ in terms of the position $\dot{x}(t)$.
 
 #### a- for  $\alpha=0, \gamma=0$ (Harmonic oscillator)
 
This is the usual harmonic oscillator, with the position and the velocity oscillating for ever, as you
can see on this plot:
 
 
 
 The kinetic energy $K$ and the potential energy $U$ both oscillate for ever, and the total energy $E = K + U$ is equal to a constant at any time. This is conservation of energy.
 
 
 
 
 The phase portrait is closed, and the trajectories go for ever around this ellipse representing pairs
of points $x(t)$ and $\dot{x}(t)$. We see that $x(t)$ and $\dot{x}(t)$ never reach $0$.





```python
x_0 =5
x_dot_0 =0
gamma = 0
w_0 = 0.8
alpha =0
phi =10
a=0
b=60
H=0.001
x,v,t=Oscillator_func(x_0,x_dot_0,gamma,w_0,alpha,phi)
K = (1/2)*((w_0)**2)*(x**2)
U = (1/2)*(v**2)
plt.figure(figsize=(8,14))
plt.subplot(3,1,1)
plt.plot(t,x)
plt.plot(t,v)
plt.legend(['x(t)','v(t)'])
plt.subplot(3,1,2)
plt.plot(t,K)
plt.plot(t,U)
plt.legend(['Kinetic Energy','Potential Energy'])
plt.subplot(3,1,3)
plt.plot(x,v)
plt.show()
```

#### b- for  $\alpha=0, \gamma=0.1$ (Damped oscillator)


In this case we have $\gamma^2 < \omega_{0}^2$, which corresponds to the under-damped regime discussed in class. The solution has an oscillatoray behavior with an exponential damping, as you can see on this
plot:


Es expected, the mechanical energy is not conserved because of friction (this is why the motion
eventually stops). You can see on this plot that energy goes to zero:



In the phase portrait the trajectory starts at $x_{0}=5$ and $\dot{x}_{0}=0$, which is our initial condition, and
is attracted to the point $x = 0$ and $\dot{x}=0$. This is consistent with what we see above on the plots
of $x(t)$ and $\dot{x}(t)$.




```python
x_0 =5
x_dot_0 =0
w_0 = 0.8
alpha =0
phi =10
gamma = 0.1
x,v,t=Oscillator_func(x_0,x_dot_0,gamma,w_0,alpha,phi)
plt.figure(figsize=(10,16))
plt.subplot(3,1,1)
plt.plot (t,x)
plt.plot(t,v)
plt.legend(['x(t)','v(t)'])
plt.subplot(3,1,2)
plt.plot(t,K)
plt.plot(t,U)
plt.legend(['Kinetic Energy','Potential Energy'])
plt.subplot(3,1,3)
plt.plot(x,v)
plt.show()
```

#### c- for  $\alpha=1, \gamma=0.1$ (Damped oscillator with driving) 

In this case the motion is initially damped, as in the previous case, but after some time the driving force takes over and becomes dominant, so the motion is again oscillatory, where the oscillations are dictated by the driving. This can be seen on the plot.


The energy is initially decreasing, as in the damped case, but then the driving starts to inject energy in the system, so energy increases once again and starts to oscillate between kinetic and potential.


In the phase portrait the trajectory starts at $x_{0}=5$ and $\dot{x}_{0}=0$, which is our initial condition. It is initially attracted to the center, because of the damping as in the previous case, but then we see that the trajectory does not converge towards $(0, 0)$, but goes back to a stable circle.



```python
x_0 =5
x_dot_0 =0
gamma = 0.1
w_0 = 0.8
alpha =1
phi =10
x,v,t=Oscillator_func(x_0,x_dot_0,gamma,w_0,alpha,phi)
K = (1/2)*(w_0)**2*x**2
U = (1/2)*v**2
plt.figure(figsize=(10,16))
plt.subplot(3,1,1)
plt.plot (t,x)
plt.plot(t,v)
plt.legend(['x(t)','v(t)'])
plt.subplot(3,1,2)
plt.plot(t,K)
plt.plot(t,U)
plt.legend(['Kinetic Energy','Potential Energy'])
plt.subplot(3,1,3)
plt.plot(x,v)
plt.show()
```


```python
{
 "cells": [],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 2
}


```




    {'cells': [], 'metadata': {}, 'nbformat': 4, 'nbformat_minor': 2}




```python

```

f7debb9dcc38d79dffb01f2246f9543f1a93f68b
ipynb
Jupyter Notebook
audrey_demafo_CM.ipynb
AudreyDemafo/AudreyDemafo.github.io
ccc2ed8898fc87d073c2ebfc4b136ccda21f4103
MIT
audrey_demafo_CM.ipynb
AudreyDemafo/AudreyDemafo.github.io
ccc2ed8898fc87d073c2ebfc4b136ccda21f4103
MIT
audrey_demafo_CM.ipynb
AudreyDemafo/AudreyDemafo.github.io
ccc2ed8898fc87d073c2ebfc4b136ccda21f4103
MIT
Qwen/Qwen-72B
1. YES
2. YES