The Harmonic oscillation under damping and driven force

The Content

• 1. The external force
• 2. No friction and driven
• 3. No driven with finite friction
  • 3.1 Over Damping $\gamma>2k$
  • 3.2 Under Damping $\gamma<2k$
  • 3.3 Critical Damping $\gamma=2k$
  • 3.4 Comparing different damping
• 4. With finite damping and driven
  • 4.1 The resonance
  • 4.2 The numerical solution

Loading the package

using Plots;
using LaTeXStrings
gr();

1. The external force

Defining the force including the elastic force, frictious force and driven force

F(x,v,t)=-k*x-γ*v+Fd*cos(ωₛ * t);

2. No friction and driven

Initial the parameters

x0=2;     #The initial position
V0=0;      #The initial velocity
m=1;      #The mass of the particle
k=2;      #The elastic coefficient
γ=0;     #The frictious force
Fd=0;     #The driven force
ωₛ=2;    #The driven frequency

Δt=0.001;t=0:Δt:6;     #The simulation time

xt=zeros(1,length(t));  #Initialization of the position vs time xt(t)
xt[1]=x0;

Vt=zeros(1,length(t));  #Initialization of the velocity vs time Vt(t)
Vt[1]=V0;

We can get the trajectory by numerical way

for i=1:length(t)-1
    a=F(xt[i],Vt[i],t[i])/m;    #The acceleration
    Vt[i+1]=Vt[i]+a*Δt;        #The velocity in next step
    xt[i+1]=xt[i]+Vt[i]*Δt+a*Δt^2/2;            #The numerical trajectory
end

Meanwhile, we can also obtain the analytic solution

xp(t)=x0*cos(sqrt(k/m)*t);        #Analytic solution

Comparing the numerical and analytic solution

plot(Array(t),xt',w=2,c=:black,lab="Eular")
plot!(xp,0,6,w=2,c=:red,lab="Exact")

3. No driven with finite friction

Initial the parameters

x0=2;     #The initial position
V0=0;      #The initial velocity
m=1;      #The mass of the particle
k=1;      #The elastic coefficient
Fd=0;     #The driven force
ωₛ=2;    #The driven frequency

3.1 Over-Damping $\gamma>2k$

Initial the friction

γ=6;     #The frictious force

Δt=0.0001;t=0:Δt:6;     #The simulation time

xt=zeros(1,length(t));  #Initialization of the position vs time xt(t)
xt[1]=x0;

Vt=zeros(1,length(t));  #Initialization of the velocity vs time Vt(t)
Vt[1]=V0;

We can get the trajectory by numerical way

for i=1:length(t)-1
    a=F(xt[i],Vt[i],t[i])/m;    #The acceleration
    Vt[i+1]=Vt[i]+a*Δt;        #The velocity in next step
    xt[i+1]=xt[i]+Vt[i]*Δt+a*Δt^2/2;            #The numerical trajectory
end

Meanwhile, we can also obtain the analytic solution

β=γ/(2*m);ω₀=√(k/m);
α₊=-β+√(β^2-ω₀^2);α₋=-β-√(β^2-ω₀^2);
A=α₋*x0/(α₋-α₊);B=α₊*x0/(α₊-α₋);

xp(t)=A*exp(α₊*t)+B*exp(α₋*t);        #Analytic solution

Comparing the numerical and analytic solution

plot(Array(t),xt',w=2,c=:black,lab="Eular")
plot!(xp,0,6,w=2,c=:red,lab="Exact")

3.2 Under-Damping $\gamma<2k$

Initial the friction

γ=0.5;     #The frictious force

Δt=0.0001;t=0:Δt:12;     #The simulation time

xt=zeros(1,length(t));  #Initialization of the position vs time xt(t)
xt[1]=x0;

Vt=zeros(1,length(t));  #Initialization of the velocity vs time Vt(t)
Vt[1]=V0;

We can get the trajectory by numerical way

for i=1:length(t)-1
    a=F(xt[i],Vt[i],t[i])/m;    #The acceleration
    Vt[i+1]=Vt[i]+a*Δt;        #The velocity in next step
    xt[i+1]=xt[i]+Vt[i]*Δt+a*Δt^2/2;            #The numerical trajectory
end

Meanwhile, we can also obtain the analytic solution

β=γ/(2*m);ω₀=√(k/m);ωp=√(ω₀^2-β^2);

xp(t)=x0*exp(-β*t)*(cos(ωp*t)+β/ωp*sin(ωp*t));        #Analytic solution

Comparing the numerical and analytic solution

plot(Array(t),xt',w=2,c=:black,lab="Eular")
plot!(xp,0,12,w=2,c=:red,lab="Exact")

3.3 Critical-Damping $\gamma=2k$

Initial the friction

γ=2*k;     #The frictious force

Δt=0.0001;t=0:Δt:12;     #The simulation time

xt=zeros(1,length(t));  #Initialization of the position vs time xt(t)
xt[1]=x0;

Vt=zeros(1,length(t));  #Initialization of the velocity vs time Vt(t)
Vt[1]=V0;

We can get the trajectory by numerical way

for i=1:length(t)-1
    a=F(xt[i],Vt[i],t[i])/m;    #The acceleration
    Vt[i+1]=Vt[i]+a*Δt;        #The velocity in next step
    xt[i+1]=xt[i]+Vt[i]*Δt+a*Δt^2/2;            #The numerical trajectory
end

Meanwhile, we can also obtain the analytic solution

β=γ/(2*m);ω₀=√(k/m);

xp(t)=x0*exp(-β*t)*(1+β*t);        #Analytic solution

Comparing the numerical and analytic solution

plot(Array(t),xt',w=2,c=:black,lab="Eular")
plot!(xp,0,12,w=2,c=:red,lab="Exact")

3.4 Comparing different damping

Defining the damping function for all case

function x_damp(γ,k,m,t,x0)
    β=γ/(2*m);ω₀=√(k/m);
    if β> ω₀   #Over damping
        α₊=-β+√(β^2-ω₀^2);α₋=-β-√(β^2-ω₀^2);
        A=α₋*x0/(α₋-α₊);B=α₊*x0/(α₊-α₋);
        x_damp=A*exp.(α₊.*t).+B*exp.(α₋.*t);
    elseif β < ω₀ #Under damping
        ωp=√(ω₀^2-β^2);
        x_damp=x0*exp.(-β.*t).*(cos.(ωp.*t)+β/ωp*sin.(ωp.*t));
    else           #Critical damping
        x_damp=x0*exp.(-β.*t).*(1 .+β.*t);
    end
end

x_damp (generic function with 1 method)

Changing the damping coefficent, we can see the fast damping is critical damping

plot(t->x_damp(0,1,1,t,1),0,12,w=2,c=:black,ls=:dash,lab="No damping")
plot!(t->x_damp(0.2,1,1,t,1),0,12,w=2,c=:red,lab="under damping \\gamma=0.2")
plot!(t->x_damp(1,1,1,t,1),0,12,w=2,c=:red,ls=:dash,lab="under damping \\gamma=1")
plot!(t->x_damp(1.5,1,1,t,1),0,12,w=2,c=:red,ls=:dot,lab="under damping \\gamma=1.5")
plot!(t->x_damp(2,1,1,t,1),0,12,w=2,c=:black,lab="critical damping \\gamma=2")
plot!(t->x_damp(2.5,1,1,t,1),0,12,w=2,c=:blue,ls=:dash,lab="over damping \\gamma=2.5")
plot!(t->x_damp(3,1,1,t,1),0,12,w=2,c=:blue,lab="over damping \\gamma=3")

4. With finite damping and driven

4.1 The resonance

The resonance function is $x(t)=|Z_0|Cos(\omega_s t-\psi)$, where $tan(\psi)=\frac{2\omega_s\beta}{\omega_0^2-\omega_s^2}$ and the amplitude $|Z_0|$ is

Z₀(ωₛ)=Fd/m/√((ω₀^2-ωₛ^2)^2+4*ωₛ^2*β^2)

Z₀ (generic function with 1 method)

Its maximum $Z_m(\beta)$ exists at $\omega_s=\sqrt{\omega_0^2-2\beta^2}$ with value

Zₘ(β)=Fd/(2*m*β*√(ω₀^2-β^2));ωₘ(β)=√(ω₀^2-2*β^2);

Initial the parameters

ω₀=1;Fd=1;m=1;

The amplitude vs $\omega_s$

β=0.1;plot(Z₀,0.,1.5,w=2,c=:red,lab=L"\beta=0.1");
β=0.15;plot!(Z₀,0,1.5,w=2,c=:blue,lab=L"\beta=0.15");
β=0.2;plot!(Z₀,0,1.5,w=2,c=:yellow,lab=L"\beta=0.2");
β=0.25;plot!(Z₀,0,1.5,w=2,c=:cyan,lab=L"\beta=0.25");
β=0.3;plot!(Z₀,0,1.5,w=2,c=:black,lab=L"\beta=0.3");
plot!(ωₘ,Zₘ,0.1,0.5,w=2,c=:black,ls=:dash,lab=L"Z_m")
plot!(xlabel=L"\omega_s",ylabel=L"Z_0")

4.2 The numerical solution

Initial the parameters

x0=2;     #The initial position
V0=0;      #The initial velocity
m=1;      #The mass of the particle
k=1;      #The elastic coefficient
γ=0.2;     #The frictious force
Fd=1;     #The driven force
ωₛ=1;    #The driven frequency

Δt=0.001;t=0:Δt:12;     #The simulation time

xt=zeros(1,length(t));  #Initialization of the position vs time xt(t)
xt[1]=x0;

Vt=zeros(1,length(t));  #Initialization of the velocity vs time Vt(t)
Vt[1]=V0;

We can get the trajectory by numerical way

for i=1:length(t)-1
    a=F(xt[i],Vt[i],t[i])/m;    #The acceleration
    Vt[i+1]=Vt[i]+a*Δt;        #The velocity in next step
    xt[i+1]=xt[i]+Vt[i]*Δt+a*Δt^2/2;            #The numerical trajectory
end

Meanwhile, we can also obtain the analytic solution

xp(t)=x0*cos(sqrt(k/m)*t);        #Analytic solution

Comparing the numerical and analytic solution

plot(Array(t),xt',w=2,c=:black,lab="Euler")
plot!(xp,0,12,w=2,c=:red,lab="No friction and driven")