Forced Vibrations : Mathematica Explorations

Today we will use Mathematica to explore some of the properties surrounding Forced Vibrations. If the differential equation is linear, we can solve it by hand using the techniques for second order linear constant coefficient differential equations.

If the differential equation is nonlinear, we will use Mathematica's built-in command NDSolve to solve the differential equation. This is a bit of a "black-box" at this point, but before the end of the semester we will see how to solve equations like this by converting them to systems and using a modification of Taylor's method of Order n.

A Note About Mathematica Syntax

I am using the following structure to define and solve the DEs since I want to be able to change the parameters m,γ, k easily.

L[y_,m_,γ_,k_]:=m D[y,{t,2}] + γ D[y,{t,1}] + k y 

Clear[u]
sol=First[u/.DSolve[{L[u[t],1,1,5/4]==3Cos[1*t],u[0]==2,u'[0]==3},u,t]];
sol[t]

Note you can, in your own work, choose not to use the L[u[t],1,1,5/4], and just enter the DE directly. Then your solution would look like

Clear[u]
sol=First[u/.DSolve[{u''[t]+u'[t]+5/4 u[t]==3Cos[1*t],u[0]==2,u'[0]==3},u,t]];
sol[t]

For numerical solutions using NDSolve, you can do something similar.

Clear[u]
sol=First[u/.NDSolve[{u''[t]+u'[t]+5/4 u[t]==3Cos[1*t],u[0]==2,u'[0]==3},u,{t,0,10}]];
sol[t]

A Linear Example: Steady - State and Transient Part of Solutions

L[y_,m_,γ_,k_]:=m D[y,{t,2}] + γ D[y,{t,1}] + k y

Let's solve an IVP based on this DE using Mathematica's DSolve command.

Clear[u]
sol=First[u/.DSolve[{L[u[t],1,1,5/4]==3Cos[1*t],u[0]==2,u'[0]==3},u,t]];
sol[t]
Plot[sol[t],{t,0,20},ImageSize->{500,Automatic},PlotRange->{All,All},PlotStyle->Black,BaseStyle->Medium,AxesLabel->{"t","u(t)"}]

The transient is the part of the solution that approaches zero as t->∞.

The part that does not go to zero is called the steady-state part of the solution.

steadystate[t_]=(1/17) (12  Cos[t] Cos[2 t]+51 Sin[t]+3  Cos[2 t] Sin[t]-3  Cos[t] Sin[2 t]+12  Sin[t] Sin[2 t]);

transient[t_]=(1/17) (E^(-t/2)) (22 Cos[t]+14 Sin[t]);

Plot[{steadystate[t],transient[t],sol[t]},{t,0,20},PlotStyle->{{Thick,Dashed},{Thick,Blue},{Thick,Black}},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"t","u(t)"}]

We could be more precise if needed, but just by eye it looks like the transient is not contributing very much to the solution by the time we get to t>10.

If we want to roughly determine the amplitude of the steady state solution, we can plot a line and adjust until we get the amplitude correct by eye.

Plot[{steadystate[t],transient[t],sol[t],2.91},{t,0,20},PlotStyle->{{Thick,Dashed},{Thick,Blue},{Thick,Black}},ImageSize->{500,Automatic},PlotRange->{All,All},AxesLabel->{"t","u(t)"}]

From the above, it looks like the amplitude of the steady state solution is roughly 2.91.

How Amplitude R of Steady State-Solution Depends on Frequency of External Driving Force

It is possible to derive an equation for how the amplitude depends on the frequency of the steady-state solution (Section 3.8 Equation (13) in 10th Edition of Boyce & DiPrima).

R=(((Subscript[F, 0])/k)/(Sqrt[(((1-((((ω/(Subscript[ω, 0]))))^2)))^2)+Γ((((ω/(Subscript[ω, 0]))))^2)]))

In our case above, the parameters are m=1,γ=1,k=5/4,(Subscript[F, 0])=3,ω=1,(Subscript[ω, 0])=(Sqrt[k/m])=(Sqrt[5])/2, Γ=(γ^2)/(m k)=4/5. Substituting these values into Eq. (13) we get

R[ω_]=(1/(Sqrt[(((1-((((ω^2)2/(Sqrt[5])))^2)))^2)+(4/5)*((((ω^2)2/(Sqrt[5])))^2)]))

Plot[R[ω],{ω,0,2},AxesLabel->{"ω","R"},BaseStyle->Medium,PlotStyle->Thick]

So we see that the maximum amplitude of the steady state solution occurs when ω=1, and then drops of fairly rapidly.

Phase Plots: (du/dt)vs u

In physics and chemistry, it is often useful to plot the position against the velocity. These types of plots are called phase plots. We can generate these in Mathematica by creating a table of values and plotting the points.

data=Table[{sol'[t],sol[t]},{t,0,20,0.02}];

Hookesphase1=ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"}]

Notice that in the phase plot, we see that the "trajectory" looks like it is approaching an ellipse, which we call a limit cycle. This is one of the common possibilities you see for a linear differential equation. The ellipse is what you approach as the transient solution approaches zero. Previously, it looked like the transient wasn't contributing much for t>10, so we can plot that region to see the steady-state behaviour of the system.

data=Table[{sol'[t],sol[t]},{t,10,20,0.02}];

Hookesphasesteadystate=ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},BaseStyle->Medium,PlotRange->{All,All},AxesLabel->{"u(t)","u'(t)"}]

A Nonlinear Example: Hardening Spring

A hardening spring does not obey Hooke's law; it is nonlinear, and the effect of a hardening spring is that the restoring force is stronger for larger displacements than for a spring that satisfies Hooke's law. For the model below, if ϵ>0 we have a hardening spring (if ϵ<0 we have a softening spring).

Clear[L,u,sol,t,y]
L[y_,m_,γ_,k_,e_]:=m D[y,{t,2}] + γ D[y,{t,1}] + k y +e*(y^3)

Because this is a fairly complicated nonlinear DE, DSolve is unable to find a solution. We must resort to a numerical solution, using the NDSolve command. The result of NDSolve is an interpolating function solution that you can work with in Mathematica as if you had an explicit function solution u(t). You can plot, take derivatives, integrate, etc. the interpolating function solution.

sol=First[u/.NDSolve[{L[u[t],1,1,5/4,0.1]==3Cos[1*t],u[0]==2,u'[0]==3},u,{t,0,100}]];
sol[t]
Plot[sol[t],{t,0,100},ImageSize->{500,Automatic},BaseStyle->Medium,PlotRange->{All,All},PlotStyle->Black,AxesLabel->{"t","u(t)"}]

Notice that the transient part of the solution does contribute early on, but the solution ultimately settles down into something that is periodic once t>20 (again, we could be more precise in this estimation if needed). Since we have an interpolating function solution, we cannot easily split it up into a transient and steady-state part.

Phase Plots

Notice the phase plot still approaches a limit cycle, but it is more complex that what we saw for the linear case.

data=Table[{sol'[t],sol[t]},{t,0,100,0.02}];

Hardeningphase1=ListPlot[data,Joined->True,PlotStyle->{Thick,Red},ImageSize->{500,Automatic},PlotRange->{All,All},AxesLabel->{"u(t)","u'(t)"}]

Here is the steady-state solution once the transient is contributing very little; notice that it is no longer an ellipse.

data=Table[{sol'[t],sol[t]},{t,20,100,0.02}];

Hardeningphasesteadystate=ListPlot[data,Joined->True,PlotStyle->{Thick,Red},ImageSize->{500,Automatic},PlotRange->{All,All},AxesLabel->{"u(t)","u'(t)"}]

Compare the steady-state solutions from the Hooke's law spring to the Hardening spring:

Show[Hookesphasesteadystate,Hardeningphasesteadystate]

This makes sense, since a hardening spring will restrict the motion of the spring more than a Hooke's law spring.

Sensitivity of Steady-State to Initial Conditions

If we want to examine how the steady-state solution for the hardening spring depends on the initial conditions, we need to do some numerical experimentation. We will choose a variety of initial conditions, and generate numerical solutions. This is a different differential equation than we worked with above (different parameters).

list={};
Table[

sol=First[u/.NDSolve[{L[u[t],1,2,3,1]==4Cos[2t]Sin[3t],u[0]==i,u'[0]==i},u,{t,0,100}]];

data=Table[{sol'[t],sol[t]},{t,0,100,0.02}];

p[i]=ListPlot[data,Joined->True,PlotStyle->{Thick,Hue[i/5]},BaseStyle->Medium,ImageSize->{500,Automatic},PlotRange->{All,All},AxesLabel->{"u(t)","u'(t)"}];

list=Append[list,p[i]];
,{i,1,5}];

Show[list,PlotRange->{All,All}]

They all end up in the same steady-state regardless of the initial condition.  This is very clear if we only plot what the solutions look like for t>20. Notice the steady-state phase space plot is fairly complicated.

list={};
Table[

sol=First[u/.NDSolve[{L[u[t],1,2,3,1]==4Cos[2t]Sin[3t],u[0]==i,u'[0]==i},u,{t,0,100}]];

data=Table[{sol'[t],sol[t]},{t,20,100,0.02}];

p[i]=ListPlot[data,Joined->True,PlotStyle->{Thick,Hue[i/5]},BaseStyle->Medium,ImageSize->{500,Automatic},PlotRange->{All,All},AxesLabel->{"u(t)","u'(t)"}];

list=Append[list,p[i]];
,{i,1,5}];

Show[list]

A Linear Example with Periodic Driving Function

Often we have a periodic forcing function, and therefore it is useful to know how to program those in Mathematica. Here's one way of getting it done. The Floor command is the greatest integer function.

Periodic Function on t>0:

g(t)=f(t), 0≤ t ≤ T; 
g(t+T) = g(t) where T >0 is a constant.

In this example, we are using f(t) =sin(t) (e^(-t/5)) and a period of T=3π.

f[t_]=Sin[t]*Exp[-t/5]
g[t_,T_]:=f[t-T*Floor[(t/T)]]
Plot[g[t,3Pi],{t,0,9π},PlotStyle->{Thick,Red},BaseStyle->Medium,AxesOrigin->{0,0},AxesLabel->{"t","u(t)"}]

Clear[L,u,sol,t,y]
L[y_,m_,γ_,k_,e_]:=m D[y,{t,2}] + γ D[y,{t,1}] + k y +e*(y^3)

sol=First[u/.NDSolve[{L[u[t],1,0.2,6.4,0]==g[t,3Pi],u[0]==0,u'[0]==0},u,{t,0,100}]];

p1=Plot[sol[t],{t,0,100},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,PlotStyle->Black,AxesLabel->{"t","u(t)"}]

Phase Plots

Here is the steady-state phase plot:

data=Table[{sol'[t],sol[t]},{t,80,100,0.02}];

ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"}]

A Linear Example with Square Wave Driving Function

Square wave function:

Clear[g,t,T]
f[t_]=Piecewise[{{6,0<t<1},{0,1<t<2}}];
g[t_,T_]:=f[t-T*Floor[(t/T)]]
Plot[g[t,2],{t,0,10},PlotStyle->{Thick,Red},BaseStyle->Medium,AxesOrigin->{0,0}]

Here's another way of doing this:

Clear[g]
g[t_,A_,T_]:=(A/2)*((-1)^Floor[t/T]+1)
Plot[g[t,6,1],{t,0,10},PlotStyle->{Thick,Red},BaseStyle->Medium,AxesOrigin->{0,0}]

Clear[L,u,sol,t,y,g,t,T]
L[y_,m_,γ_,k_,e_]:=m D[y,{t,2}] + γ D[y,{t,1}] + k y +e*(((y))^3)

f[t_]=Piecewise[{{6,0<t<1},{0,1<t<2}}];
g[t_,T_]:=f[t-T*Floor[(t/T)]]

sol=First[u/.NDSolve[{L[u[t],1,0.2,0.2,0]==g[t,2],u[0]==0,u'[0]==0},u,{t,0,100}]];
sol[t]
p1=Plot[sol[t],{t,0,100},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,PlotStyle->Black,AxesLabel->{"t","u(t)"}]

Phase Plots: u' vs u

data=Table[{sol'[t],sol[t]},{t,0,100,0.02}];

ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"}]

data=Table[{sol'[t],sol[t]},{t,90,100,0.02}];

ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"}]

Problem: Nonlinear

Solve the IVP and sketch u(t) as well as the phase plane plot and the steady state phase space plot.

12 u'' + 3u' +2u+(u^5)=R(t), u(0)=1, u'(0)=-2,
where 
 R(t)=20cos(2t)sin(t), 0 ≤ t ≤ π, is periodic with period π.

Barry's Solution

Define the periodic forcing function :

Clear[g,t,T]
f[t_]=20Cos[2t]Sin[t]
g[t_,T_]:=f[t-T*Floor[(t/T)]]
Plot[g[t,Pi],{t,0,24},PlotStyle->{Thick,Red},BaseStyle->Medium,AxesOrigin->{0,0}]

Solve the differential equation:

sol=First[u/.NDSolve[{12 u''[t] + 3u'[t] +2u[t]+((u[t])^5)==g[t,Pi],u[0]==1,u'[0]==-2},u,{t,0,100}]];
sol[t]
p1=Plot[sol[t],{t,0,100},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,PlotStyle->Black,AxesLabel->{"t","u(t)"}]

Phase space plot:

data=Table[{sol'[t],sol[t]},{t,0,100,0.02}];

ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"}]

Steady state phase space plot:

data=Table[{sol'[t],sol[t]},{t,60,100,0.02}];

ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"},AxesOrigin->{0,0}]

Problem: Linear with Resonance

Solve the IVP and sketch u(t) as well as the phase plane plot. This one can be solved by hand.

u''+u=(1/2)cos(t), u(0)=1, u'(0)=-2,

Barry's Solution

Solve the differential equation (note this can be solved for u(t)):

sol=First[u/.DSolve[{u''[t]+u[t]==1/2Cos[t],u[0]==1,u'[0]==-2},u,t]];
sol[t]
p1=Plot[sol[t],{t,0,100},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,PlotStyle->Black,AxesLabel->{"t","u(t)"}]

Phase space plot:

data=Table[{sol'[t],sol[t]},{t,0,100,0.02}];

ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"}]

The steady-state solution does not approach a closed curve. This is an example of resonance, where the frequency of the driving force matches the frequency of the oscillator.

Problem: Linear with Beats

Solve the IVP and sketch u(t) as well as the phase plane plot. This one can be solved by hand.

u''+u=(1/2)cos(0.9t), u(0)=1, u'(0)=-2,

Barry's Solution

Solve the differential equation (note this can be solved for u(t)):

sol=First[u/.DSolve[{u''[t]+u[t]==1/2*Cos[0.9t],u[0]==1,u'[0]==-2},u,t]];
sol[t]
p1=Plot[sol[t],{t,0,300},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,PlotStyle->Black,AxesLabel->{"t","u(t)"}]

Phase space plot:

data=Table[{sol'[t],sol[t]},{t,0,300,0.02}];

ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"}]

The steady-state solution is actually quite complicated here. This is an example of beats.

data=Table[{sol'[t],sol[t]},{t,100,300,0.02}];

ListPlot[data,Joined->True,PlotStyle->{Thick,Blue},ImageSize->{500,Automatic},PlotRange->{All,All},BaseStyle->Medium,AxesLabel->{"u(t)","u'(t)"}]