This material includes the derivation of harmonic oscillator solutions using Sage in the following classic cases:
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$.
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.
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:
load('cas_utils.sage')
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:
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:
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:
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
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
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)$:
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$).
#@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))
free_oscilator_xt(w0=1, x0=1)
free_oscilator_xt(w0=3, x0=1)
free_oscilator_xt(w0=1/2, x0=1)
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:
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:
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:
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.
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}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)
In Sage the above solution contains two constants which have no "handles":
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:
[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:
_K1
Since the solution is in the form of a periodic function multiplied by exponent, we might want to draw the envelope:
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:
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.
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)) }
showmath(phi_osc.subs(sub3))
Alternatively, you can match the solution to a bit more complicated pattern:
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)) }
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:
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:
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))
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:
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)
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:
var('t')
plot( phi_anal.subs({_K1:-1,_K2:1,omega0:1,g:2.1}), (t,0,15), figsize=(10,2))
var('t')
plot( phi_anal.subs({_K1:-1,_K2:0,omega0:1,g:2.1}), (t,0,15), figsize=(10,2))
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:
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.
var('t')
plot( phi_anal.subs({_K1:-1,_K2:1,g:1}), (t,0,15), figsize=(10,2))
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$:
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}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:
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.
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.
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.
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:
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)
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$.
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:
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):
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:
var('a b x')
assume(a>0)
(a*sin(x)+b*cos(x)).subs(sub3).show()
and expand the obtained formula to check:
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:
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:
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.
expr = r_sz.numerator().collect(sin(omega*t))
showmath(expr)
expr = expr.subs(sub3)
showmath(expr)
expr = expr.subs(sub3a)
showmath(expr)
and in all its glory our formula is presented as:
r_szczegolne = (expr/expr_denom).canonicalize_radical()
showmath(r_szczegolne)
If for some reason we would like to pick its coefficients then we have:
w0 = SR.wild(0)
w1 = SR.wild(1)
m = r_szczegolne.match(w0*sin(w1))
showmath(m)
check:
showmath(w0.subs(m)*sin(w1.subs(m)))
Thus, the amplitude of the oscillator in the limit $t\to\infty$ is:
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.
phase = w1.subs(m)-omega*t-pi
showmath(phase)
The maximal amplitude depends on the attenuation coefficient:
#@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))
plot_A_phase(g_=0.1),plot_A_phase(g_=1)
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:
var('phid',latex_name=r'\dot\varphi')
var('phi',latex_name=r'\varphi')
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_}
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')
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()
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$:
showmath(A.diff(omega))
sol = solve(A.diff(omega), omega)
showmath(sol)
omega_max = sol[1].rhs()
showmath(omega_max)
showmath( (omega_max^2).factor() )
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.
pars = {a:1,omega0:1}
plot3d( A.subs(pars), (omega,0,3),(g,0.2,1) ).show(viewer='tachyon',figsize=4)
pars = {a:1,omega0:1}
contour_plot( A.subs(pars), (omega,0,3),(g,0.2,1) )
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:
showmath( r_szczegolne )
Dissipated power is equal to: force times the velocity, so:
P = g*(r_szczegolne.diff(t))^2
We want to integrate the power over one period of oscillations
showmath( P.integrate(t,0,2*pi) )
We can plot this formula together with an expression for amplitude.
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