# Harmonic oscillator with computer algebra

## Harmonic oscillator
This material includes the derivation of harmonic oscillator solutions using Sage in the following classic cases:

- free oscillator
- free damped oscillator
- forced oscillator with damping: resonance phenomenon
 

### What is a harmonic oscillator?


The harmonic oscillator is a point particle with mass $m$ in the force field, which depends linearly on the position. The restoring force to the equilibrium position depends linearly on the amount of the deflection. An example of such a force may be the reaction of a flexible body applying Hook's law. In other words, we have a point particle of  mass $m$ in a square potential: 

$$ U(x) = \frac{1}{2}k x^2$$
The  Newton equation for this system reads:
$$ m a = -\frac{dU(x)}{dx} = - kx,$$
and because acceleration $a$ is the second derivative of the position with respect to time one gets the differential equation:
$$ m \ddot x   =  -kx.$$
This is the equation of motion for the harmonic oscillator. Another method of obtaining it is to substitute the force derived from the linear span - that is the Hook's law $\vec F = -k \vec x$, to equate $m\vec a = \vec F$.

### Small oscillations

The harmonic oscillator is an extremely important model that appears in many real situations, both in classical and quantum systems. We now show why? Note that for any potential $U(x)$, which at some point $x_0$, say $x_0=0$,  has a local minimum, its Taylor series takes the form:

$$ U(x) = U(0) + U'(0)x + \frac{1}{2}U''(0) x^2+ ...$$

Since we have assumed that at $x_0=0$ the potential has a local minimum, the derivative at this point is zero, i.e.,  $U'(0) = 0$ and therefore:

$$ U(x) = U(0) + \frac{1}{2}U''(0) x^2+ ...$$

The Taylor's series can be approximated by keeping only the first two non-vanishing terms. The first term $U(0)$ is constant and can be neglected because the corresponding force does not depend on this constant.  As a result, the potential is approximated by the expression:

$$ U(x) = \frac{1}{2}U''(0) x^2.$$

The potential is identical to the potential of the harmonic oscillator and  the second derivative of the potential at the point of its minimum $U"(0)$ can be identified with the elastic constant $k$. So we see that the movement of the particle in the vicinity of the potential minimum  can be approximated by a harmonic oscillator. This fact is the reason why the harmonic oscillator is such frequently used in physics.



## Free oscillator

Consider the Newton equation for a harmonic oscillator in a vacuum (there is no dissipation,no a frictional force) It has a form:
$$  m \ddot x   =  -kx.$$
It is conveniet to convert this equation to the form 
$$  \ddot x   =  -\omega_0^2x,$$
where $\omega_0 = \sqrt{\frac{k}{m}}$ is a positive number.

It is a linear second order  ordinary differential equation with constant coefficients. Equations of this class can be easily solved - assuming the form of a solution and substituting it into the Newton  equation. With the computer algebra system included in SageMath, we can simplify this procedure using the desolve function:


```python
load('cas_utils.sage')
```


```python
var('omega0 x0')
assume(omega0>0)

var('t')
X = function('X')(t)

osc = diff(X,t,2) == -omega0^2*X 
showmath( osc )
```

Because we use a set of variables for operations in expressions containing derivatives, the above cell gives us a mathematical formula of the differential equation representing the harmonic oscillator.

With this expression, we can use computer algebra to solve the differential equation:


```python
phi_anal = desolve(osc,dvar=X,ivar=t,show_method=True)
showmath(phi_anal)
```

First, we see that we must assume that $\omega_0$ is non-zero. Otherwise, the solution would have a different form.

Secondly, we see that the so-called a general solution depends on two constants $K_1$ and $K_2$. We can determine these two constants knowing the position and velocity of the oscillator at some given time $t_0$. Because the equation of motion do not depend explicitly on time, the solution $x(t) = x(t, t_0) = x(t-t_0) depends on the time difference $t-t_0. Therefore, without the limitation of generality, $t=0$ can be used as an initial moment. 


Sage can also determine constants. If at the moment $t=0$  the oscillator is  in the rest state: $x(0)=x_0$ and $v(0)=0$, then we have:


```python
phi_anal = desolve(osc,dvar=X,ivar=t,ics=[0,x0,0])
showmath(phi_anal)
```

If, at $t=0$, the initial conditions are:  $x(0)=0$ and $v(0)=v_0$, then we have:


```python
var('v0 x0')
phi_anal = desolve(osc,dvar=X,ivar=t,ics=[0,0,v0])
showmath(phi_anal)
```

Finally, for a general case when  $x(0)=x_0$ and $v(0)=v_0$ the result reads


```python
var('v0 x0')
phi_anal = desolve(osc,dvar=X,ivar=t,ics=[0,x0,v0])
showmath(phi_anal)
```

We can also assume specific numerical values $x(0)=2$ and $v(0)=3$. Then 


```python
phi_anal = desolve(osc,dvar=X,ivar=t,ics=[0,2,3])
showmath(phi_anal)
```

Let's plot this solution with initial conditions $x_0=2$  $v_0=3$ for  $\omega_0=1$. We see that the solution is tangent to $3/\omega_0$ line that passes the point $(0,2)$:


```python
plot(phi_anal.subs(omega0==1),(t,0,5),figsize=4)+\
 plot( 3*t+2,(t,-1,1),color='gray')+\
 point([0,2],color='red')
```

We can examine how the solution depends on the parameters $\omega_0$ and $x_0$ (with $v_0=0$).


```python
#@interact
def free_oscilator_xt(w0=slider(0.1,3,0.1,default=1),x0=slider(0.1,3,0.1,default=1)):
    phi_anal = desolve(osc,dvar=X,ivar=t,ics=[0,1,0])
    p = plot(phi_anal.subs({omega0:1}),(t,0,7),color='gray')
    phi_anal = desolve(osc,dvar=X,ivar=t,ics=[0,x0,0])
    p += plot(phi_anal.subs({omega0:w0}),(t,0,7))
    p.show( ticks=[[0,pi/2,pi,3/2*pi,2*pi],1/2],tick_formatter="latex",figsize=(6,2))
```


```python
free_oscilator_xt(w0=1, x0=1)
free_oscilator_xt(w0=3, x0=1)
free_oscilator_xt(w0=1/2, x0=1)
```

## Damped oscillator
As we have seen from previous considerations, the free harmonic oscillator exhibits oscillations  infinitively long. In real situations,  systems are in contact with its surroundings (environment) which introduces a dissiption process. In other words, the system is subjected to damping forces, which cause the oscillations  to decay in time. It is the so called damped oscillator. A special case of the friction force, which is taken into account when testing the motion of the oscillator - is the frictional force  proportional to the velocity:
$$ \vec F_D = - \Gamma \vec v, $$
where $\Gamma$ is a damping (friction) constant. It is so called the Stokes force or the drag force
Such friction occurs, for example, when the ball is in the liquid with small Reynolds numbers and is called viscous friction. If the Reynolds number is large, then the friction depends in a more complicated way on the velocity. E.g. in aerodynamics,  a quadratic relationship is a good model.

Consideration of a linearly rate-dependent friction force has the basic advantage that, despite adding an additional term  to the equation of motion, we still  deal with a linear differential equation with constant coefficients and a full mathematical analysis of solutions can be carried out. So we have the equation:


$$m \ddot x   + \Gamma \dot x  + kx  =  0.$$

Let's divide this equation by $m$ and introduce the rescaled damping constant $\gamma=\Gamma/m$  and 
 $\omega_0 =k/m$. Then 

$$\ddot x   + \gamma \dot x  + \omega_0^2 x  =  0.$$

From  theory of linear differential equations we know that we should consider the roots of a characteristic polynomial  for a given equation and depending on their nature we obtain various different solutions. It can be easily done assuming that the solution of the Newton equation in an exponent $e^{k t}$ and substitute this form into ODE:


```python
var('k t')
var('omega omega0')
var('g', latex_name='\gamma')

f = exp(k*t)
eq = f.diff(t,2)+g*f.diff(t)+omega0^2*f
showmath(eq)
```

Now we can factorize this expression:


```python
eq.factor().show()
```

It is zero for all $t$ only when ${\gamma} k + k^{2} + \omega_{0}^{2}=0$. This is a quadratic equation which has the following solutions:


```python
showmath( (eq.factor()*exp(-k*t)).solve(k) )
```

There are three different values of the determinant:

$$\Delta = \sqrt{\gamma^2-4 \omega_0^2}$$

for which the solution has qualitatively different properties:

 - Two complex roots - damped oscillations.
 - Two real negative roots - strongly damped motion, no oscillation.
 - One degenerate root  ($\Delta=0$ ) - critical vibrations and  it can have one maximum.

### Two complex roots:
Since both $\gamma$ and $\omega_0$ are positive and the expression for $\Delta$  factorizes:

$$ \gamma^2-4 \omega_0^2 = ( \gamma-2 \omega_0)( \gamma+2 \omega_0),$$ 

the expression under square root will be negative if and only if 

\begin{equation}
\label{eq:osc_1}
\gamma-2 \omega_0 < 0
\end{equation}



```python
var('omega omega0')
var('g', latex_name='\gamma')

forget()
assume(g-2*omega0<0)
assume(g>0)
assume(omega0>0)
show( assumptions() ) 

osc = diff(X,t,2) == -g*diff(X,t)-omega0^2*X 

showmath(osc)
```


```python
phi_anal = desolve(osc,dvar=X,ivar=t)
showmath(phi_anal)
```

In Sage the above solution contains two  constants which have no "handles":


```python
phi_anal.variables()
```

It means that expession like `.subs({_K1:2})` will throw an exception.

We can however define symbolic variables for all constants in solution. Note that those constants start with  underscore character, thus we can do it automatically:


```python
[var(str(s)) for s in phi_anal.variables() if str(s).startswith("_")]
```

Now we have `_K1` variable available, and we can use it in substitutions:


```python
_K1
```

Since the solution is in the form of a periodic function multiplied by exponent, we might want to draw the envelope:


```python
var('t')
pars={_K1:11,_K2:2,omega0:1,g:.51}
A = sqrt(_K1^2+_K2^2).subs(pars)
plot( phi_anal.subs(pars), (t,0,15), figsize=(10,2)) + \
    plot( (A*exp(-1/2*g*t)).subs(pars), (t,0,15),color='gray' ) +\
    plot( (-A*exp(-1/2*g*t)).subs(pars), (t,0,15),color='gray' )
```

**Where did the formula for A come from?**

We can use the trigonometric identity  which allow to add $\sin$ and $\cos$ functions with the same frequency: 
$$\sqrt{2} \sin\left(\frac{1}{4} \, \pi + {x}\right)=\cos\left({x}\right) + \sin\left({x}\right)$$


$$\sqrt{a^{2} + b^{2}} \sin\left({x} + \arctan\left(\frac{b}{a}\right)\right)$$

We have:



```python
phi_osc = phi_anal.coefficient(e^(-1/2*g*t))
showmath(phi_osc)
```

We can transform the linear combination of the $\sin$ and $\cos$ functions into one function with a different amplitude. In order to accomplish this in Sage, we can use  wildcard substitutions.


```python
w0 = SR.wild(0)
w1 = SR.wild(1)
w2 = SR.wild(2)
sub3 =  { w1*sin(w0)+w2*cos(w0):sqrt(w1^2+w2^2)*sin(w0+arctan(w2/w1)) }
```


```python
showmath(phi_osc.subs(sub3))
```

Alternatively, you can match the solution to a bit more complicated pattern:


```python
w0 = SR.wild(0)
w1 = SR.wild(1)
w2 = SR.wild(2)
w3 = SR.wild(3)
sub4 =  { w3*(w1*sin(w0)+w2*cos(w0)):w3*sqrt(w1^2+w2^2)*sin(w0+arctan(w2/w1)) }
```


```python
showmath(phi_anal.subs(sub4))
```

On the other hand if have expression in the form of  $\sin(x+\phi) e^{-x}$, then the exponent term is an envelope: 


```python
plot(sin(x)*exp(-x/3) ,(x,0,10),figsize=(6,2))+\
 plot([exp(-x/3),-exp(-x/3)],(x,0,10),color='red')
```

we have therefore: 


```python
a,b = 2,1.2
plot([(a*sin(x)+b*cos(x))*exp(-x/3),sqrt(a^2+b^2)*exp(-x/3)],\
     (x,0,10), figsize=(6,2))
```

### Two real roots

If the determinant $\Delta$ is positive (but non-zero), i.e.:

\begin{equation}
\label{eq:osc_2}
\gamma-2 \omega_0 > 0,
\end{equation}

then we have: 


```python
var('omega omega0')
var('g', latex_name='\gamma')
forget()
assume(g-2*omega0>0)
assume(g>0)
assume(omega0>0)
show( assumptions() ) 
osc = diff(X,t,2) == -g*diff(X,t)-omega0^2*X 
showmath(osc)
```


```python
phi_anal = desolve(osc,dvar=X,ivar=t)
showmath(phi_anal)
```

We see that the solution does not contain periodic functions but  is the sum of two exponents with negative arguments.

Task: Prove that the argument of the exponent 
$k_{1} e^{\left(-\frac{1}{2} \, {\left(\gamma - \sqrt{\gamma^{2} - 4 \, \omega_{0}^{2}}\right)} t\right)}$,
for $t>0$, is negative.

In this case, solutions are decaying without oscillations, however, for certain parameters it is possible that a single extremum will occur:



```python
var('t')
plot( phi_anal.subs({_K1:-1,_K2:1,omega0:1,g:2.1}), (t,0,15), figsize=(10,2))
```


```python
var('t')
plot( phi_anal.subs({_K1:-1,_K2:0,omega0:1,g:2.1}), (t,0,15), figsize=(10,2))
```

### The degenerate case

Consider the case when the determinant is zero, 

\begin{equation}
\label{eq:osc_3}
\gamma-2 \omega_0 = 0
\end{equation}

Let us have a look how the general solution looks:


```python
var('omega omega0')
var('g', latex_name='\gamma')

forget()
assume(g-2*omega0==0)
assume(g>0)
assume(omega0>0)
show( assumptions() ) 
osc = diff(X,t,2) == -g*diff(X,t)-omega0^2*X 
showmath(osc)
phi_anal = desolve(osc, dvar=X, ivar=t)
showmath(phi_anal)
```

This case is called:  critically damped. 


```python
var('t')
plot( phi_anal.subs({_K1:-1,_K2:1,g:1}), (t,0,15), figsize=(10,2))
```

### Dissipated power

If we know analytical solutions then  one can evaluate  how much of system energy is dissipated to the environment. For a given value of  $\omega_0$ we might expect that power remains constant in two cases: $\gamma=0$ and $\gamma\to\infty$. Let's check it out - first we can obtain analytical solutions in oscilating and damped regimes. We will take initial condition $x_0=1$ and $v_0=1$:


```python
forget()
assume(g>2*omega0)
assume(g>0)
assume(omega0>0)
x_damped = desolve(osc, dvar=X, ivar=t, ics=[0,1,0] )

forget()
assume(g<2*omega0)
assume(g>0)
assume(omega0>0)
x_oscil = desolve(osc, dvar=X, ivar=t, ics=[0,1,0] )
```

Now, let's define functions which characterize dissipated power $P$  (it is energy or work per unit time):

\begin{equation}
\label{eq:power_harm}
P = F_{friction} v = (\gamma v) v =  \gamma v^2
\end{equation}

and total energy (for $m=1$):


\begin{equation}
\label{eq:E_harm}
E_{tot} = \frac{1}{2}v^2 + \frac{1}{2} \omega_0 x^2 
\end{equation}




```python
E = lambda x: 1/2 * x.diff(t)^2 + 1/2*omega0*x^2
P = lambda x: g * (x.diff(t))^2
```

Now we can plot power  and trajectory as a function of time:


```python
plt_power = plot([P(x_damped).subs({omega0:1,g:g_}) for g_ in [2.01,4,10]],\
     (t,0,10), color='black')
plt_power += plot([P(x_oscil).subs({omega0:1,g:g_}) for g_ in [.1,1,1.99]],\
      (t,0,10), color='red')
plt_power.show(figsize=(6,2))
```

In the above figure, the dissipated power  $P$ is shown for selected values of the damping constant. 


```python
plt_E = plot([E(x_damped).subs({omega0:1,g:g_}) for g_ in [2.01,4,10,1e3]],\
             (t,0,10),color='black')
plt_E += plot([E(x_oscil).subs({omega0:1,g:g_}) for g_ in [1e-3,.1,1,1.99]],\
              (t,0,10),color='red')
plt_E.show(figsize=(6,2))
```

In the above figure, the total energy $E_{tot}$ is depicted for various damping constants. 


```python
plt_x = plot([(x_damped).subs({omega0:1,g:g_}) for g_ in [2.01,4,10]],\
             (t,0,10),color='black')
plt_x += plot([(x_oscil).subs({omega0:1,g:g_}) for g_ in [.1,1,1.99]],\
              (t,0,10),color='red')
plt_x.show(figsize=(6,2))
```

The above figure shows the solutions $x(t)$ for different damping constant. 

## Forced harmonic oscillator
Now, we take into account the driving in the form of a time-dependent periodic force $\sin(\omega t)$.   In this case, the equation of motion  contains a term independent of the position $x(t)$, but explicitely dependent  on time. Such an equation is a linear non-homogeneous differential equation and we can give its analytical  solutions.
Let's see how they look:


```python
var('a omega omega0')
var('g', latex_name='\gamma')
forget()
assume(g-2*omega0<0)

Phi = function('Phi')(t)
assume(g>0)
assume(omega0>0)
osc = diff(Phi,t,2)+ g*diff(Phi,t) + omega0^2*Phi -a*sin(omega*t)

showmath(osc)
```


```python
phi_anal,method = desolve(osc,dvar=Phi,ivar=t,show_method=True)
print(method)
showmath(phi_anal)
```

It can be seen that we have a solution which is a sum of two terms: a general solution of a free damped oscillator, hence a homogeneous equation (contains $k_1$ and $k_2$ constants) and a solution that does not contain free constants. It is  called a special solution to the heterogeneous equation. It can be noted that this solution is the one that survives in the limit $t\to\infty$.


```python
r_sz = phi_anal.operands()[1]
showmath(r_sz)
```

Let's transform a special solution to the  form: $$A\sin(\omega t+\phi).$$ 
To do this, we extract  the numerator and denominator:


```python
expr_denom = r_sz.denominator()
expr = r_sz.numerator()
showmath([expr,expr_denom])
```

The numerator is a sum of $\sin$ and $\cos$ with different amplitudes. The formula: $$a \sin\left(x\right) + b \cos\left(x\right) =\sqrt{a^{2} + b^{2}} \sin\left(x + \arctan\left(\frac{b}{a}\right)\right)$$ can be used
for this purpose to put with wildcards (an. Wildcards):


```python
w0 = SR.wild(0)
w1 = SR.wild(1)
w2 = SR.wild(2)
w3 = SR.wild(3)
sub3 =  { w1*sin(w0)+w2*cos(w0):sqrt(w1^2+w2^2)*sin(w0+arctan(w2/w1)) }
```

Let's check how the substitution works on the formula:


```python
var('a b x')
assume(a>0)
(a*sin(x)+b*cos(x)).subs(sub3).show()
```

and expand the obtained formula to check:


```python
assume(a>0)
(a*sin(x)+b*cos(x)).subs(sub3).full_simplify().show()
```

Everything looks nice, but we have a classic problem to choose the function's branch. Let's look at the graph of the $\tan(x)$ function:


```python
print (arctan(10)-pi).n(),(arctan(-10)-pi).n()

plot(tan(x),(x,-4.,4),detect_poles='show',figsize=5,ymin=-10,ymax=10).show()
sum([plot(arctan(x)+n*pi,(x,-10.,10),figsize=(5,3)) for n in range(-1,2)])
```

We should take the branch $y\in -\pi .. 0$ and the function in Sage takes $y\in -\pi ..\pi$. Therefore, it is better to use the two-argument `arctan2`:


```python
w4 = SR.wild(4)
w5 = SR.wild(5)
sub3a = {arctan(w4/w5):(arctan2(w4,w5)-pi)}
```

Returning to our oscillator, we can use the substitution:

*Note* `collect (sin (omega * t))` is needed for the pattern to be in a form in which the substitution can be used.


```python
expr = r_sz.numerator().collect(sin(omega*t))
showmath(expr)
```


```python
expr = expr.subs(sub3)
showmath(expr)
```


```python
expr = expr.subs(sub3a)
showmath(expr)
```

and in all its glory our formula is presented as:


```python
r_szczegolne = (expr/expr_denom).canonicalize_radical()
showmath(r_szczegolne)
```

If for some reason we would like to pick its coefficients then we have:


```python
w0 = SR.wild(0)
w1 = SR.wild(1)
m = r_szczegolne.match(w0*sin(w1))
showmath(m)
```

check:


```python
showmath(w0.subs(m)*sin(w1.subs(m)))
```

Thus, the amplitude of the oscillator in the limit $t\to\infty$  is:


```python
A = -w0.subs(m)
showmath(A)
```

It is a famous formula related to the phenomenon of resonance. If the damping is weak then in the case of $\omega_0 \to \omega$, the amplitude of the special solution will increase to very large values, and will be infinite if the system is not damped. This means that the energy from the external source can be pumped in particularly efficient way  if the frequency of the external force is consistent with the eigen-frequency of the system.

The strength of this phenomenon was seen by soldiers passing over the bridge in Angers (http://en.wikipedia.org/wiki/Angers_Bridge), which had just its own frequency equal to the frequency of the military march. The bridge during the march of the column of the army, stared to swing to such an amplitude that it was destroyed.


```python
phase = w1.subs(m)-omega*t-pi
showmath(phase)
```

The maximal amplitude depends on the attenuation coefficient:


```python
#@interact
def plot_A_phase(g_=slider(0.01,1,0.01,label="$\gamma$",default=0.2)):
    pars = {g:g_,a:1,omega0:1}
    print ( g^2*omega^2 + omega^4 - 2*omega^2*omega0^2 + omega0^4 ).subs(pars)
    p = plot(A.subs(pars),(omega,0,10))
    p += plot(phase.subs(pars),(omega,0,10),color='red',detect_poles="show")
    p.show(figsize=(6,2))
```


```python
plot_A_phase(g_=0.1),plot_A_phase(g_=1)
```

### Numerical analysis

We can also numerically integrate the equation of motion for the harmonic oscillator. It requires to know all parameters and initial conditions numerically before, and it does not allow for easy analysis of the generic properties of the solutions. On the other hand numerical procedure can be applied easily to system of ODE which do not have analytical solution. Here, let us compare results obtained from numerical solution with $t\to\infty$ formula:


```python
var('phid',latex_name=r'\dot\varphi')
var('phi',latex_name=r'\varphi')

```


```python
rhs = solve(osc,diff(Phi(t), t, 2))[0].rhs()
#@interact
def plot_A_traj(g_=slider(0.01,2.2,0.01,label="$\gamma$",default=0.2),\
                omega_=slider(0.01,5.134,0.01,label="$\omega$",default=1.4)):
    pars = {g:g_,a:1,omega0:1}
    print ( g^2*omega^2 + omega^4 - 2*omega^2*omega0^2 + omega0^4 ).subs(pars)
    p = plot(A.subs(pars),(omega,0,10),figsize=(10,2))
    p += plot(phase.subs(pars),(omega,0,10),color='red',detect_poles="show")
    pars = {omega:omega_,omega0:1,a:1,g:g_}
    
```


```python
pars = {omega:1+0.1,omega0:1,a:0.2,g:.2}
ode = [phid,rhs.subs({Phi:phi,diff(Phi,t):phid}).subs(pars)]
times = srange(0,80,0.001)
ics = [0.0,2.1]
sol = desolve_odeint(ode,ics,times,[phi,phid])
line( zip(times,sol[::1,0]),figsize=(10,2) ) + \
 plot( r_szczegolne.subs(pars), (t,0,80),color='red') +\
 plot( (a*sin(omega*t)).subs(pars), (t,0,80),color='gray')
```


```python
rhs = solve(osc,diff(Phi(t), t, 2))[0].rhs()
#@interact
def plot_A_traj(g_=slider(0.01,2.2,0.01,label="$\gamma$",default=0.2),\
                omega_=slider(0.01,5.134,0.01,label="$\omega$",default=1.4)):
    pars = {g:g_,a:1,omega0:1}
    print ( g^2*omega^2 + omega^4 - 2*omega^2*omega0^2 + omega0^4 ).subs(pars)
    p = plot(A.subs(pars),(omega,0,10),figsize=(10,2))
    p += plot(phase.subs(pars),(omega,0,10),color='red',detect_poles="show")
    pars = {omega:omega_,omega0:1,a:1,g:g_}
    ode = [phid,rhs.subs({Phi:phi,diff(Phi,t):phid}).subs(pars)]
     
    times=srange(0,80,0.001)
    ics=[0.0,2.1]
    sol = desolve_odeint(ode,ics,times,[phi,phid])
    r_szczegolne2 = a*sin(omega*t +pi- \
                          arctan2(-g*omega,(omega^2 - omega0^2)))/\
           sqrt(g^2*omega^2+ omega^4 - 2*omega^2*omega0^2 + omega0^4)
    p2 = line( zip(times,sol[::1,0]),figsize=(10,2) ) + \
     plot( r_szczegolne.subs(pars), (t,0,80),color='red') +\
     plot( (a*sin(omega*t)).subs(pars), (t,0,80),color='gray')
    p += point([omega_,0])
    p.show()
    p2.show()
```


```python
plot_A_traj(g_=0.2,omega_=1.)
plot_A_traj(g_=4.2,omega_=1.)
```

We can further analyze the formula for resonance and e.g. find its maximum in $omega$:


```python
showmath(A.diff(omega))
```


```python
sol = solve(A.diff(omega), omega)
showmath(sol)
```


```python
omega_max = sol[1].rhs()
showmath(omega_max)
```


```python
showmath( (omega_max^2).factor() )
```


```python
p=[]
for g_ in srange(0.1,1,0.1)+srange(1,1.41,.2):
    pars = {g:g_,a:1,omega0:1}
    omega_max_v =  omega_max.subs(pars).n()
    if omega_max_v.is_real():
        p.append( point( (omega_max_v,A.subs(pars).subs(omega==omega_max_v) ),color='red' ) ) 
    p.append(  plot(A.subs(pars),(omega,0,3)) ) 
sum(p).show(ymax=11)
```

In above formula we see that in the weak damping limit the resonance frequency tends to internal frequency of the oscillator. But with increased damping we position moves towards lower values of frequency. Also it can be easily seen that for 

$$ -\gamma^2+2 \omega_0^2 <0$$ 

because:

$$ -\gamma^2+2 \omega_0^2 = (\sqrt 2 \omega_0-\gamma)(\sqrt 2 \omega_0+\gamma)$$ 

for 

$$ \sqrt 2 \omega_0 < \gamma, $$

the resonance disappears - i.e. there is not maxima in $A(\omega)$ dependence. 


```python
pars = {a:1,omega0:1}
plot3d( A.subs(pars), (omega,0,3),(g,0.2,1) ).show(viewer='tachyon',figsize=4)
```


```python
pars = {a:1,omega0:1}
contour_plot( A.subs(pars), (omega,0,3),(g,0.2,1) )
```

###  Absorbed power 


We can also analyze how much power is absorbed by the periodically driven harmonic oscillator. We can use the solution $t\to\infty$ which has a form:


```python
showmath( r_szczegolne )
```

Dissipated power is equal to: force times the velocity, so:


```python
P = g*(r_szczegolne.diff(t))^2
```

We want to integrate the power over one period of oscillations


```python
showmath( P.integrate(t,0,2*pi) )
```

We can plot this formula together with an expression for amplitude. 


```python
p=[]
for g_ in [0.25,0.5,1,2,4]:
    pars = {g:g_,a:1,omega0:1}
    assume(omega>0)
    omega_max_v =  omega_max.subs(pars).n()
    if omega_max_v.is_real():
        p.append( point( (omega_max_v,A.subs(pars).subs(omega==omega_max_v) ),color='red' ) ) 
    p.append(  plot(A.subs(pars),(omega,0,3), color='gray') ) 
    p.append(  plot(P.subs(pars).integrate(t,0,2*pi/omega)/2,(omega,0,3),color='red') ) 
sum(p).show(figsize=(6,3))
```

\newpage

0c3784f20378aef10baa8210cb03d8d1dd23dd6f
ipynb
Jupyter Notebook
01-Harmonic_oscillator.ipynb
marcinofulus/Mechanics_with_SageMath
6d13cb2e83cd4be063c9cfef6ce536564a25cf57
MIT
01-Harmonic_oscillator.ipynb
marcinofulus/Mechanics_with_SageMath
6d13cb2e83cd4be063c9cfef6ce536564a25cf57
MIT
2022-01-30T16:45:58.000Z
2022-01-30T16:45:58.000Z
01-Harmonic_oscillator.ipynb
marcinofulus/Mechanics_with_SageMath
6d13cb2e83cd4be063c9cfef6ce536564a25cf57
MIT
2020-11-15T08:26:14.000Z
2022-02-12T13:07:16.000Z
Qwen/Qwen-72B
1. YES
2. YES