Chapter 7: Oscillations

7.1 Simple Harmonic Oscillator

Consider a block of mass m attached to a spring on a horizontally frictionless table. Expriments have shown that the force exerted by the spring is linearly proportional to the distance x the spring is displaced from its equilibrium position. Mathematically, this is

F=-k x

where k is the spring constant. Its value depends on the particular spring used.
The set of first-order ODEs that describes this situation for the position x(t) and velocity u(t) of the spring is therefore

m(du/dt)=-k x	(1)

(dx/dt)=u		(2)

with initial conditions

x(0)=(Subscript[x, i])
u(0)=(Subscript[u, i])

Note thas this is a coupled set of ODEs because equation 1 canot be solved independently of Equation 2 since it depends on x. In previous cases examined with falling bodies and projectiles we could solve for the velocity independently of the position. This is not the case here.
A general solution to this problem can be found analytically. It is given by:

Remove["Global`*"]

(* "E" is for Exact solution *)
{xE[t_],uE[t_]}={x[t],u[t]}/.Simplify[DSolve[{m u'[t]==-k x[t],x'[t]==u[t],u[0]==ui,x[0]==xi},{x[t],u[t]},t]]//First

Since there is no dissipation, i.e. no friction, the total mechanical energy of this system should be conserved. Let's check this. The total mechanical energy is the sum of the potential U and kinetic K energy. The potential energy is given by:
	U=(1/2)(kx^2)
and the kinetic energy is
	K=(1/2)(mu^2)
Adding these contributions, we should get a constant.

EnergyE=(1/2)k ((xE[t])^2)+(1/2)m ((uE[t])^2)//Simplify

Of course, this problem can also be solved numerically as long as the values of the parameters and boundary conditions are defined. Here's an algorithm using the midpoint scheme and the following parameters: 
	(Subscript[x, i])=0.1 m
k=0.2 kg/(s^2)
(Subscript[u, i])=0
m=0.5 kg

For the purpose of this example we'll use a timestep Δt=0.1 s and the simulation period will be 50 s. Here's sample code to achieve this along with a visual comparison with the exact solution found above.

(* Parameters *)
xi=0.1; k=0.2; ui=0;m=0.5;
T=50;  (* Simulation period *)
Δt=0.1; (* time step*)

(* Initial condition *)
t[1]=0;
u[1]=ui;
x[1]=xi;

nmax=Round[T/Δt];

(* Midpoint scheme for u and x *)
Do[t[n+1]=n Δt;
sx=u[n];
su=-((k x[n])/m);
xtmp=((Δt)/2)sx+x[n];  (* Note that since u depends on x we need to update x here *)
utmp=((Δt)/2)su+u[n];
sx=utmp;
su=-((k xtmp)/m);
x[n+1]=Δt sx+x[n];
u[n+1]=Δt su+u[n],
{n,1,nmax}]

xN=Interpolation[Table[{t[n],x[n]},{n,1,nmax+1}]]; (* xN for Numerical sol. of x *)
uN=Interpolation[Table[{t[n],u[n]},{n,1,nmax+1}]]; (* uN for Numerical sol. of u *)

Plot[{xN[t],xE[t]},{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]
Plot[{uN[t],uE[t]},{t,0,T},AxesLabel->{"t(s)","u(m/s)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]

It would be illuminating to examine how the error between the numerical algorithm and the exact solution evolves in time. Let simply examime the difference between the two solutions. To do this we will create a list like this:

Plot[xE[t]-xN[t],{t,0,T},AxesLabel->{"t(s)","Error (m)"}]
Plot[uE[t]-uN[t],{t,0,T},AxesLabel->{"t(s)","Error (m/s)"}]

This shows that the error grows in time. The same for the velocity u.

Let's plot the total energy E for the numerical solution to see how constant it is:

(* Total energy from the numerical solution *)
EnergyN[t_]:=(1/2)k ((xN[t])^2)+(1/2)m ((uN[t])^2)

(* Comparative graphs *)
Plot[{EnergyE,EnergyN[t]},{t,0,T},AxesLabel->{"t(s)","Energy (J)"},PlotStyle->{{RGBColor[0,0,1],Thickness[0.01]},{RGBColor[1,0,0],Thickness[0.01]}}]

As expected, the total energy for the exact solution is constant and conserved. However, the total energy of the numerical solution is not conserved and increases with time. This is due to the errors in the numerical algorithm. While the graph above may suggest that the energy is increasing a lot it would be more illuminating to examine how the energy increase relative to the initial energy of the system. Let's instead examine the percentage change in the energy over the course of the simulation

Plot[100 Abs[((EnergyN[t]-EnergyN[0])/(EnergyN[0]))],{t,0,T},AxesLabel->{"t(s)","% Energy change"},PlotStyle->{RGBColor[1,0,0],Thickness[0.02]}]

So, the maximum energy change is really only 0.2% over this simulation. For a practical application this may be acceptable.

In this example the quantity k/m has the units of (s^(-2)) which is units of angular frequency squared. Therefore, the equation for harmonic equation is often written more generally as

(((d^2)x)/(dt^2))+(Subscript[ω, 0]^2)x=0		(Equation of simple harmonic oscillator)

where (Subscript[ω, 0])=(Sqrt[k/m]) is the characteristic angular frequency.  The frequency (in Hz, or cycles per second) is f=((Subscript[ω, 0])/(2 π)) This is a second-order differential equation that, as before, can be reduced to a set of 2 first-order ODEs, i.e.

(du/dt)+(Subscript[ω, 0]^2)x=0		
	(dx/dt)=u

Note: the form of the solution given by Mathematica is  B Cos[ω t] + C Sin[ω t]

sol=DSolve[{u'[t]==-(ω^2) x[t],x'[t]==u[t],x[0]==xi,u[0]==ui},{x[t],u[t]},t]

X[t_]=x[t]/.sol[[1]]//Simplify

This is equivalent to the form A Cos[ω t + ϕ], since

A Cos[ω t+ϕ]//TrigExpand

and with A=(Sqrt[(xi^2)+((((ui/ω)))^2)])  and tan(ϕ)=-(ui/(xi ω))
The parameter A is the amplitude and ϕ is the phase angle.  In terms of the amplitude, the energy is simply (1/2)k (A^2)=(1/2)m (ω^2) (A^2).  Let’s check this.

A=(Sqrt[(xi^2)+((((ui/ω)))^2)])

(1/2)m (ω^2)  A^2//Simplify

(1/2)m X'[t]^2+(1/2)m (ω^2) X[t]^2//Simplify

7.2 Damped Harmonic Oscillator

The linearly damped oscillator is described with

(((d^2)x)/(dt^2))+γ(dx/dt)+(Subscript[ω, 0]^2)x=0		(Equation for damped oscillations)

The 2nd order ODE above can be reduced to a system of 2 first-order ODEs, i.e.,

(du/dt)+γu+(Subscript[ω, 0]^2)x=0		
	(dx/dt)=u

Where γ is a damping coefficient that depends on the type of oscillator. With initial conditions, i.e.,

x(0)=(Subscript[x, i])
u(0)=(Subscript[u, i])

this system has a general analytic solution

Remove["Global`*"]

{xE[t_],uE[t_]}={x[t],u[t]}/.Simplify[DSolve[{u'[t]+γ u[t]+((ω0)^2)x[t]==0,x'[t]==u[t],x[0]==xi,u[0]==ui},{x[t],u[t]},t]]//First

Case 1: Underdamping

Where are the Sines and Cosines in the above expression for the general solution?  They are hidden in the exponential.  When (γ^2)<4((Subscript[ω, 0])^2), (Sqrt[(γ^2)-4((Subscript[ω, 0])^2)]) becomes imaginary and, as you recall, Exp[i θ] = Cos[θ]+ i Sin[θ].  In this case, the solution reduces to

xUnder[t_]:=A Exp[-(γ/2) t]Cos[(Subscript[ω, 1]) t + ϕ]

where

(Subscript[ω, 1])=ω0(Sqrt[1-((γ^2)/(4 ((ω0)^2)))]);
ϕ=ArcTan[-(ui/(xi (Subscript[ω, 1])))-(γ/(2 (Subscript[ω, 1])))];
A=(Sqrt[(xi^2)+((((2 ui+xi γ))^2)/(4 ((Subscript[ω, 1])^2)))]);

Note that this damped oscillator has a smaller frequency ((Subscript[ω, 1])< (Subscript[ω, 0])), and hence a longer period.

What about the energy of the oscillator in this case?

Energy[t_]=(1/2)k ((xUnder[t])^2)+(1/2)m ((xUnder'[t])^2)

This looks a little complicated, but in the case of light damping, i.e. when γ/4 << (Subscript[ω, 0])  and hence (Subscript[ω, 1])~(Subscript[ω, 0])=(Sqrt[k/m]), this simplifies approximately to

EnergyApprox[t_]:=(1/2)k (A^2) Exp[-γ t]

Let’s confirm all of this by considering the following parameters.

(* Parameters *)
xi=0.1; k=0.2; ui=2;m=0.5;
γ=0.1; (* Damping factor *)
ω0=(Sqrt[(k/m)]); (* Frequency *)
(* Simulation period *)
T=50;

the solution of the damped harmonic oscillator looks like:

xE[t]

or with

(Subscript[ω, 1]):=ω0(Sqrt[1-((γ^2)/(4 ((ω0)^2)))]);
ϕ:=ArcTan[-(ui/(xi (Subscript[ω, 1])))-(γ/(2 (Subscript[ω, 1])))];
A:=(Sqrt[(xi^2)+((((2 ui+xi γ))^2)/(4 ((Subscript[ω, 1])^2)))]);

xUnder[t]

Plot[{xE[t],xUnder[t],A Exp[-γ t/2]},{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{{Red,Thick},Black,Green}]

And the energies

Plot[{Energy[t],EnergyApprox[t]},{t,0,T},AxesLabel->{"t(s)","Energy(J)"},PlotStyle->{{Red,Thick},Black}]

We see that the energy decays exponentially.  The characteristic decay time τ is the time required for the oscillator’s energy to fall to (e^(-1))  of its initial value.  Since (approximately)

	E[t]=(1/2)k (A^2) Exp[- γ t], then we have
	E[τ]=(1/2) k (A^2) Exp[-1]=(1/2) k (A^2) Exp[-τ γ] and hence τ = 1/γ

Overdamping

As the damping parameter γ increases, the motion of the oscillator dies down more quickly...to a point.

γ=0.05
Plot[{xUnder[t],A Exp[-γ t/2]},{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{Black,Green}]

γ=0.5
Plot[{xUnder[t],A Exp[-γ t/2]},{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{Black,Green},PlotRange->All]

γ=5;
Plot[{xUnder[t],A Exp[-γ t/2]},{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{Black,Green},PlotRange->All]

Where is the second curve A Exp[-γ t/2]?  With these parameters, A has gone imaginary and so the second curve can not be plotted.

A

What has happened to the oscillations?  In this case (γ^2)>4((Subscript[ω, 0])^2), and (Sqrt[(γ^2)-4((Subscript[ω, 0])^2)]) is real.  In this case the general solution reduces to the form

xOver[t_]:=C1 Exp[-(γ-(Sqrt[(γ^2)-4(Subscript[ω, 0])^2]))t/2]+C2 Exp[-(γ+(Sqrt[(γ^2)-4(Subscript[ω, 0])^2]))t/2]

and the oscillator’s motion decays exponentially with no oscillations.  One term decays more quickly, while the other decays more slowly.

Critical Damping

What happens γ = 2 ω0?

ω0

γ=2 ω0

xE[t]

We see that the general solution does not cover this case.  Let’s solve the harmonic oscillator equation again, this time explicitly putting in γ=2(Subscript[ω, 0]).

Remove["Global`*"]

γ=2ω0

{xE[t_],uE[t_]}={x[t],u[t]}/.Simplify[DSolve[{u'[t]+γ u[t]+((ω0)^2)x[t]==0,x'[t]==u[t],x[0]==xi,u[0]==ui},{x[t],u[t]},t]]//First;

xE[t]

We see a solution of a qualitatively different form than we’ve seen before,

	A Exp[-γ t/2] + B t Exp[-γ t/2]

You will likely see the method called “variation of parameters” that allows one to find such a solution in a class on ODE’s. For now, we will merely mention that in this case, the oscillator returns back to its equilibrium position faster than in the overdamped case.

(* Parameters *)
xi=0.1; k=0.2; ui=2;m=0.5;
ω0=(Sqrt[(k/m)]); (* Frequency *)
(* Simulation period *)
T=50;

γ

Plot[xE[t],{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{Red},PlotRange->All]

7.3 Forced and Damped Harmonic Oscillator

An oscillator can be forced. In this case the equation is

(((d^2)x)/(dt^2))+γ(dx/dt)+(Subscript[ω, 0]^2)x=F		(Forced and damped oscillations)

where F is a forcing function. For example it could be

F=(Subscript[F, 0])sin((Subscript[ω, f]) t+ϕ)

where (Subscript[ω, f]) is the forcing frequency, ϕ is the phase angle and (Subscript[F, 0]) the forcing amplitude.

A general analytical solution exists to this problem

Remove["Global`*"]

{xE[t_],uE[t_]}={x[t],u[t]}/.Simplify[DSolve[{u'[t]+γ u[t]+((ω0)^2)x[t]==F0 Sin[ωf t+ϕ],x'[t]==u[t],x[0]==xi,u[0]==ui},{x[t],u[t]},t]]//First

But it is quite complicated!

Let solve this problem numerically. Here's an algorithm using the midpoint scheme and the following parameters: 
	(Subscript[x, i])=0.0 m
γ=0.01 (s^(-1))
(Subscript[u, i])=0
(Subscript[ω, 0])=0.1 (s^(-1))
(Subscript[ω, f])=0.2 (s^(-1))
ϕ=0

For the purpose of this example We'll set Δt=0.1 s and the simulation period will be 1000 s. Here's a code to achieve this along with a visual comparison with the exact solution found above.

(* Parameters *)
xi=0.0; γ=0.01; ui=0;k=0.015;m=1.0;ωf=0.2;F0=0.2;ϕ=0°;
ω0=(Sqrt[(k/m)]);
T=1000;  (* Simulation period *)
Δt=0.01;

(* Initial condition *)
t[1]=0;
u[1]=ui;
x[1]=xi;

nmax=Round[T/Δt];

(* Midpoint scheme for u and x *)
Do[t[n+1]=n Δt;
sx=u[n];
su=F0 Sin[ωf t[n]+ϕ]-γ u[n]-((ω0)^2) x[n];
xtmp=((Δt)/2)sx+x[n];
utmp=((Δt)/2)su+u[n];
sx=utmp;
su=F0 Sin[ωf (t[n]+t[n+1])/2+ϕ]-γ utmp-((ω0)^2) xtmp;
x[n+1]=Δt sx+x[n];
u[n+1]=Δt su+u[n],
{n,1,nmax}]

xN=Interpolation[Table[{t[n],x[n]},{n,1,nmax+1}]];
uN=Interpolation[Table[{t[n],u[n]},{n,1,nmax+1}]];

Plot[{xN[t],xE[t]},{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]
Plot[{uN[t],uE[t]},{t,0,T},AxesLabel->{"t(s)","u(m/s)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]

We see that after a while, the system settles into a steady state.  Let’s compare the solution with the forcing term on the right hand side of the equation.

Plot[{xE[t],F0 Sin[ωf t+ϕ]},{t,0,1000}]

We see that in this case, the response of the oscillator is larget compare to the forcing term, but has the same frequency in the steady state.  To more clearly see this, let’s replot the functions, this time including a factor to make the comparison more apparent and focussing on later times.

Plot[{xE[t],30F0 Sin[ωf t+ϕ]},{t,600,1000}]

Additionally, we notice that the oscillator and the driving force appear to be completely out of phase (one is going up while the other is going down).  These observations motivates the form of the analytic solution presented in the next section.

Now let's look at the total energy of the system, which we determine from both the exact and the numerical solutions.

(* Total energy from the exact solution *)
EnergyE[t_]=(1/2)k ((xE[t])^2)+(1/2)m ((uE[t])^2);

(* Total energy from the numerical solution *)
EnergyN[t_]=(1/2)k ((xN[t])^2)+(1/2)m ((uN[t])^2);

(* Comparative graphs *)
Plot[{EnergyE[t],EnergyN[t]},{t,0,T},AxesLabel->{"t(s)","Energy (J)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]

Because the system is forced the energy oscillates once it reaches a steady state.  The two solutions are comparable but of course there is an error, as always, with the numerical solution.

As with the motion itself, the energy settles into a steady state.  What is the average steady state energy?  The time average of a function over some time interval is

	(Subscript[f, ave])=(1/((Subscript[t, 2])-(Subscript[t, 1])))(∫_(Subscript[t, 1])^(Subscript[t, 2]))f[t]dt.

Choosing a reasonable period for our oscillator, we determine the average steady-state energy to be

fave=(1/(1000-800))(∫_800^1000)EnergyN[t]\[DifferentialD]t//N

7.4 Resonance

Consider again the case of a forced, underdamped oscillator driven by an oscillating force (ω is the angular frequency of the driver),

(((d^2)x)/(dt^2))+γ(dx/dt)+(Subscript[ω, 0]^2)x=(Subscript[F, 0])cos(ω t)

for which we propose a solution of the form

x(t)= A cos(ω t + ϕ) + B (e^(-(γ/2)t))cos((Subscript[ω, 1])t+ψ)	(7.4.1)

You will recognize the rightmost term as the solution to the damped oscillator (satisfying  (((d^2)x)/(dt^2))+γ(dx/dt)+(Subscript[ω, 0]^2)x=0), with (Subscript[ω, 1])=(Subscript[ω, 0])(Sqrt[1-((γ^2)/(4 ((Subscript[ω, 0])^2)))]).  It includes the constants B and ψ that can be used to satisfy initial conditions.  Because of the exponentially decaying prefactor, this term will die out with time and leave only the term on the left as the steady-state solution, which oscillates with the same frequency as the driving force.

With some algebra, we arrive at the following expressions for the steady-state amplitude A and phase factor ϕ

A= ((Subscript[F, 0])/m)(1/(Sqrt[(((((Subscript[ω, 0])^2)-(ω^2)))^2)+(((ω γ))^2)])),
ϕ=(tan^(-1))(((γ ω)/((ω^2)-((Subscript[ω, 0])^2)))),  with ϕ between -π and 0.

We see that the amplitude of the steady state oscillations depend on the driving frequency, reaching a maximum near ω=(Subscript[ω, 0]).  Let us check these results.

Remove["Global`*"]

XE[t_]=x[t]/.DSolve[{x''[t]+γ x'[t]+((ω0)^2)x[t]==F0 Cos[ω t],x[0]==xi,x'[0]==ui},x[t],t]//First;

Our steay-state solution.

ss[t_]:=A Cos[ω t+ϕ]

A=(F0/m)(1/(Sqrt[(((((ω0)^2)-(ω^2)))^2)+(((ω γ))^2)]));
ϕ=If[ω<ω0,ArcTan[((γ ω)/((ω^2)-((ω0)^2)))],ArcTan[((γ ω)/((ω^2)-((ω0)^2)))]-π];

(The π is subtracted here to choose the correct arctan result.)

Comparing the results for some random parameters for light damping

γ=0.01; xi=1;ui=-2;k=0.25;m=1.0;ω0=(Sqrt[k/m]);ω=0.2;F0=0.2;  (* Some parameters for light damping *)

Plot[{XE[t],ss[t]},{t,0,1000}]

ω0

ω

Plot[ϕ,{ω,0.0,2 ω0},PlotRange->All]

From the graph of ϕ, we see that for a driving frequency less than ~ω0 (as we have here), the oscillator will be in phase with the driver, while for large driving frequencies, the oscillator will be out of phase.

Plot[A,{ω,0,2 ω0},PlotRange->All]

The behaviour of the steady state amplitude as a function of the driving frequency shows a large response near the natural frequency of the oscillator.  This is known as resonance.

For larger damping, the results for A and ϕ are a little different.  The response is not as sharp.

γ=0.2

Plot[ϕ,{ω,0.0,2 ω0},PlotRange->All]

Plot[A,{ω,0,2 ω0},PlotRange->All]

What is the steady-state energy of the oscillator as a funciton of the driving frequency? (Recall (Sqrt[k/m])=ω0, or k=m ((ω0)^2))

Clear[ω0,ω,γ,F0,ϕ,m,k]

ss[t]

Energy[t_]=(1/2)m ((ω0)^2) ((ss[t])^2)+(1/2)m ((ss'[t])^2)//Simplify

To get the average, we can average over 1 cycle (period is 2π/ω):

Eave=(ω/(2 π))(∫_0^(2π/ω))Energy[t]\[DifferentialD]t

In the light damping case (γ << 4(Subscript[ω, 0])), this approximates to

(Subscript[E, ave])= (1/4)(((Subscript[F, 0])^2)/m)((((ω^2)+((Subscript[ω, 0])^2)))/((((((Subscript[ω, 0])^2)-(ω^2)))^2)+(((ω γ))^2)))

(Subscript[E, ave])≈ (1/8)(((Subscript[F, 0])^2)/m)(1/((((Subscript[ω, 0])^2)-(ω^2))+(((γ/2))^2))) for light damping γ << 4(Subscript[ω, 0])

This last functional form is a Lorentzian, which is defined usually as

L[x_,a_]:=(1/π)((a/2)/((((x-x0))^2)+(((a/2))^2)))

So that it is normalized to 1 when integrated over all x

Assuming[{a>0,x0 ∈ Reals},(∫_(-∞)^∞)L[x,a]\[DifferentialD]x]

It has a maximum at x=x0, with value 2/(π a)

Maximize[L[x,a],x]

The width of the Lorentzian is a (width at half maximum)

x0=5;
Manipulate[Plot[L[x,a],{x,0,10},PlotRange-> {All,{0,10/π}}],{a,0.2,7}]

We see that as the width narrows, the maximum height increases.  As the width increases, the maximum height decreases.  Comparing to the Lorentzian our expression for the steady-state energy of the driven osciallator under light damping,

(Subscript[E, ave])≈ (1/8)(((Subscript[F, 0])^2)/m)(1/((((Subscript[ω, 0])^2)-(ω^2))+(((γ/2))^2)))

we identify the width as Δω = γ, and the maximum height as 
	(1/8)(((Subscript[F, 0])^2)/m)(4/(γ^2))
That is, the height grows as (1/((Δω)^2)).

Δω tells us how broad a range of driving frequencies is able to induce a resonant response in the oscillator.  However, we now see that the broader this range of frequency, the weaker the magnitude of the response.

Moreover, recall from section 7.2 that the characteristic decay time for an underdamped oscillator is τ = 1/γ.  Thus

		Δω τ = 1,
		
implying that if we build an oscillator the is very frequency selective (Δω small), its oscillations will decay very slowly (τ large) once we stop driving it.  On the other hand, building a broadly responsive system (Δω large) will result in a quick decay (τ small) after driving stops.

7.5 Duffing Oscillator

The general form of the equations that describe the Duffing oscillator is

(((d^2)x)/(dt^2))+γ(dx/dt)+(Subscript[ω, 0]^2)x+((βx)^3)=F		(Duffing oscillator)

where F is a forcing term, as defined in section 7.3.  This oscillator can exhibit chaotic behaviour. This problem cannot be solved analytically but can only be solved numerically. As always, the 2nd order ODE can be redued to a set of 2 first order ODEs. See the codes below that solve this problem using the midpoint scheme.

Consider the following situation (no forcing): 
	(Subscript[x, i])=1.0 m
γ=0.01 
(Subscript[u, i])=0
(Subscript[ω, 0])=1.0 (s^(-1))
(Subscript[ω, f])=(Subscript[ω, 0])/4(s^(-1))
(Subscript[F, 0])=0
ϕ=0
β=1.0 (m^(-1))(s^(-2))

Remove["Global`*"]

(* Parameters *)
xi=1.0; γ=0.01; ui=0;ω0=1.0;ωf=ω0/4;F0=0.0;β=1;
T=100;  (* Simulation period *)
Δt=0.01; (* Time step *)

(* Initial condition *)
t[1]=0;
u[1]=ui;
x[1]=xi;

nmax=Round[T/Δt];

(* Midpoint scheme for u and x *)
Do[t[n+1]=n Δt;
sx=u[n];
su=F0 Sin[ωf t[n]]-γ u[n]-((ω0)^2) x[n]-β ((x[n])^3);
xtmp=((Δt)/2)sx+x[n];
utmp=((Δt)/2)su+u[n];
sx=utmp;
su=F0 Sin[ωf (t[n]+t[n+1])/2]-γ utmp-((ω0)^2) xtmp-β (xtmp^3);
x[n+1]=Δt sx+x[n];
u[n+1]=Δt su+u[n],
{n,1,nmax}]

xN=Interpolation[Table[{t[n],x[n]},{n,1,nmax+1}]];
uN=Interpolation[Table[{t[n],u[n]},{n,1,nmax+1}]];
parametric=Interpolation[Table[{x[n],u[n]},{n,1,nmax+1}]];

Plot[xN[t],{t,0,T},AxesLabel->{"t(s)","x(m)"}]
Plot[uN[t],{t,0,T},AxesLabel->{"t(s)","u(m/s)"}]
ParametricPlot[{xN[t],uN[t]},{t,0,T},AxesLabel->{"x(m)","u(m/s)"}]

Now consider with the following parameters:

Remove["Global`*"]

(* Parameters *)
xi=1.0; γ=0.01; ui=0;ω0=1.0;ωf=ω0/4;F0=5.0;β=1;
T=200;  (* Simulation period *)
Δt=0.01; (* Time step *)

(* Initial condition *)
t[1]=0;
u[1]=ui;
x[1]=xi;

nmax=Round[T/Δt];

(* Midpoint scheme for u and x *)
Do[t[n+1]=n Δt;
s1x=u[n];
s1u=F0 Sin[ωf t[n]]-γ u[n]-((ω0)^2) x[n]-β ((x[n])^3);
xtmp=((Δt)/2)s1x+x[n];
utmp=((Δt)/2)s1u+u[n];
s2x=utmp;
s2u=F0 Sin[ωf (t[n]+t[n+1])/2]-γ utmp-((ω0)^2) xtmp-β (xtmp^3);
x[n+1]=Δt s2x+x[n];
u[n+1]=Δt s2u+u[n],
{n,1,nmax}]

xN=Interpolation[Table[{t[n],x[n]},{n,1,nmax+1}]];
uN=Interpolation[Table[{t[n],u[n]},{n,1,nmax+1}]];

Plot[xN[t],{t,0,T},AxesLabel->{"t(s)","x(m)"}]
Plot[uN[t],{t,0,T},AxesLabel->{"t(s)","u(m/s)"}]
ParametricPlot[{xN[t],uN[t]},{t,0,T},AxesLabel->{"x(m)","u(m/s)"}]

This system begins to develops a chaotic behaviour, that is a behaviour that does not repeat itself in a repeated pattern.

7.6 Pendulum

The general equation for the unforced and undamped pendulum is

(((d^2)θ)/(dt^2))+(g/L) sin(θ)=0		(Undamped and unforced pendulum)

where L is the length of the pensulum and θ the displacement angle. This second-order ODE can be replaced by the following two first-order ODEs

((Subscript[du, θ])/dt)+(g/L) sin(θ)=0
	(Subscript[u, θ])=((dθ)/dt)

where (Subscript[u, θ]) is the angular velocity of the pendulum.

This is a nonlinear set of ODEs and has a non trivial analytic solution. Let's try to solve it analytically with some initial condition, θ(0)=(Subscript[θ, i]) and (Subscript[u, θ])(0)=(Subscript[u, i])

Remove["Global`*"]

DSolve[{uθ'[t]+(g/L) Sin[θ[t]]==0,uθ[t]==θ'[t],θ[0]==θi,uθ[0]==ui},{uθ[t],θ[t]},t]

Hummm, might be a case for a numerical solution. Note however that if the pendulum has small angle the term sin{θ} ≈θ.  This can be seen by a Taylor series expansion. for example

Series[Sin[θ],{θ,0,6}]

Therefore for small angle the equation reduces to

((Subscript[du, θ])/dt)+(g/L)θ=0
	(Subscript[u, θ])=((dθ)/dt)

and has a much nicer analytical solution

DSolve[{uθ'[t]+(g/L)θ[t]==0,uθ[t]==θ'[t],θ[0]==θi,uθ[0]==ui},{uθ[t],θ[t]},t]

It will be left as an exercice, laboratory, assignment, or test question to solve numerically the pendulum problem.

7.7 Damped Pendulum

Of course the pendulum can be damped. In this case the equation is

(((d^2)θ)/(dt^2))+γ((dθ)/dt)+(g/L) sin(θ)=0		(Damped and unforced pendulum)

It will be left as an exercice, assignment or exam question to solve numerically this problem.

7.8 Damped and Forced Pendulum

Of course the pendulum can be damped and forced. In this case the equation is

(((d^2)θ)/(dt^2))+γ((dθ)/dt)+(g/L) sin(θ)=F		(Damped and forced pendulum)

where F is is some driving function.

It will be left as an exercise, assignment or exam question to solve numerically this problem.

7.9 Writing reusable routines

So far, whenever we've implemented the midpoint scheme or 4th order RK, we've taken a template and rewritten the parts that contain the right-hand-sides of our ODE's.  Wouldn't it be nice if we could use the same code over and over again?

Remove["Global`*"]

Let's start with an underdamped harmonic oscillator,

γ=1/2;
ω0=1;
T=15;

representing it as two first order ODE's.

fx[u_]:=u;
fu[x_,u_]:=-γ u - ((ω0)^2) x;

Putting in some initial conditions, we can find the analytical solution for the position and velocity.

{xE[t_],uE[t_]}={x[t],u[t]}/.Simplify[DSolve[{u'[t]==fu[x[t],u[t]],x'[t]==fx[u[t]],x[0]==1,u[0]==0},{x[t],u[t]},t]]//First

As we saw in class, the solution has exponentials with complex arguments, and so sometimes Mathematica will give a complex solution with zero as the imaginary part.

Plot[xE[t],{t,0,T},AxesLabel->{"t(s)","x(m)"}]

Now let's use NDSolve to find the solution numerically, reusing the function we used above.

{xND[t_],uND[t_]}={x[t],u[t]}/.NDSolve[{u'[t]==fu[x[t],u[t]],x'[t]==fx[u[t]],x[0]==1,u[0]==0},{x[t],u[t]},{t,0,T}]//First

Plot[xND[t],{t,0,T},AxesLabel->{"t(s)","x(m)"}]

By comparing it to the analytical solution, we see that with the default settings, NDSolve gives a solution that is accurate to about (10^(-7)) in absoulute terms.

Plot[xND[t]-xE[t],{t,0,T},AxesLabel->{"t(s)","Error"}]

Now let's write our own numerical schemes to solve the problem.  We will start with 2nd order RK (Midpoint Rule).

Δt=0.01;

It is advantageous to define as separate functions the right-hand sides of our set of first order ODE's.  This way, our RK solver is more general; all we must do is change the functions if we want to solve a different problem.

fx[u_,x_,t_]:=u;
fu[u_,x_,t_]:=-γ u - ((ω0)^2) x;

For this problem, there is no t dependence for the problem.  In general there could be and so we're leaving the door open for that.

(* Initial condition *)
t[1]=0;
u[1]=0;
x[1]=1;

nmax=Round[T/Δt];

(* Midpoint scheme for u and x *)
Do[t[n+1]=n Δt+t[1];
sx=fx[u[n],x[n],t[n]];
su= fu[u[n],x[n],t[n]];
xtmp=((Δt)/2)sx+x[n];
utmp=((Δt)/2)su+u[n];
sx=fx[utmp,xtmp,t[n]+Δt/2];
su=fu[utmp,xtmp,t[n]+Δt/2];
x[n+1]=Δt sx+x[n];
u[n+1]=Δt su+u[n],
{n,1,nmax}]

xN=Interpolation[Table[{t[n],x[n]},{n,1,nmax+1}]];
uN=Interpolation[Table[{t[n],u[n]},{n,1,nmax+1}]];

Plot[{xN[t],xE[t]},{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]
Plot[{uN[t],uE[t]},{t,0,T},AxesLabel->{"t(s)","u(m/s)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]
Plot[xN[t]-xE[t],{t,0,T},AxesLabel->{"t(s)","Error"}]

Now onto 4th order Runge-Kutta with the same time step, and using the same functions defing our set of ODE's.

(* Initial condition *)
t[1]=0;
u[1]=0;
x[1]=1;

nmax=Round[T/Δt];

(* 4th order Runge-Kutta scheme for u and x *)
Do[t[n+1]=n Δt+t[1];
s1x=fx[u[n],x[n],t[n]];
s1u= fu[u[n],x[n],t[n]];
xtmp=((Δt)/2)s1x+x[n];
utmp=((Δt)/2)s1u+u[n];
s2x=fx[utmp,xtmp,t[n]+Δt/2];
s2u=fu[utmp,xtmp,t[n]+Δt/2];
xtmp=((Δt)/2)s2x+x[n];
utmp=((Δt)/2)s2u+u[n];
s3x=fx[utmp,xtmp,t[n]+Δt/2];
s3u=fu[utmp,xtmp,t[n]+Δt/2];
xtmp=Δt s3x+x[n];
utmp=Δt s3u+u[n];
s4x=fx[utmp,xtmp,t[n+1]];
s4u=fu[utmp,xtmp,t[n+1]];

x[n+1]=((Δt)/6) (s1x + 2 s2x + 2 s3x +s4x)+x[n];
u[n+1]=((Δt)/6) (s1u + 2 s2u + 2 s3u + s4u)+u[n],

{n,1,nmax}]

xNR=Interpolation[Table[{t[n],x[n]},{n,1,nmax+1}]];
uNR=Interpolation[Table[{t[n],u[n]},{n,1,nmax+1}]];

Plot[{xNR[t],xE[t]},{t,0,T},AxesLabel->{"t(s)","x(m)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]
Plot[{uNR[t],uE[t]},{t,0,T},AxesLabel->{"t(s)","u(m/s)"},PlotStyle->{RGBColor[0,0,1],RGBColor[1,0,0]}]
Plot[xNR[t]-xE[t],{t,0,T},AxesLabel->{"t(s)","Error"}]

The above two routines can be used for other problems.  Just make sure that Δt and T are defined and the the functions fu and fx capture the right physics.

7.10 Nonlinear oscillator - perturbation approach

Remove["Global`*"]
Abig=1.2;
Asmall=0.2;
λgiven=0.3;

In this section we review the work done in class within Mathematica for working out a perturbative solution for the nonlinear oscilator we considered.  Newton's equation is given by x'' + ((Subscript[ω, 0])^2)x = λ (x^3), or
 
			x'' + ((Subscript[ω, 0])^2)x - λ (x^3)=0 

with initial conditions x(0)=A and x'(0)=0.

Let's begin by recasting the equation as  f(x)=0, where

f[x_]:= D[x,{t,2}]+((ω0)^2) x - λ (x^3)

Try a solution of the form

x=A Cos[w t]

The goal is to find w that approximately satisfies Newton's equation.  Plug this expression into f

f[x]

We need to expand (Cos^3)(w t) is terms of Cos((Subscript[n,  ])w t).  The appropriate function to use is TrigReduce, e.g.

TrigReduce[((Cos[θ])^3)]

TrigReduce[f[x]]

Now we want to collect terms that depend on Cos((Subscript[w,  ])t)

Collect[TrigReduce[f[x]],Cos[w t]]

We can make f(x)=0 approximately by finding the w that will make the coefficient of Cos(wt) zero.  We can pick parts of the expression the same way as taking elements of a list (note that the order of terms may very from what your version of Mathematica gives you).

tm=Collect[TrigReduce[f[x]],Cos[w t]][[1]]

tmp=tm[[2]]

sol=Solve[tmp==0,w]

wA=w/.sol[[2]]

This is our new frequency. It depends on the amplitude.   So our solution is A cos( w t), with w=wA as we just found.

xA= A Cos[wA t]

Let's now compare our perturbative solution to a numerical one, putting in some values for (Subscript[ω, 0]), A and λ

ω0=1;
λ=λgiven;
A= Asmall;
T=20;
x=.;
xN[t_]=y[t]/.NDSolve[{y''[t]+((ω0)^2) y[t]-λ ((y[t])^3)==0, y[0]==A,y'[0]==0},y[t],{t,0,T},PrecisionGoal->∞][[1]];
ω0
wA
Plot[{xN[t],xA,A Cos[ω0 t]},{t,0,T},PlotStyle->{{Black,Thickness[0.01]},{Red,Thickness[0.008],Dashed},{Green,Thickness[0.004]}}]

Looks like a good match.  Even the harmonic approximation is not too bad (green curve).  Let's repeat it for larger amplitude.

A=Abig;
xN[t_]=y[t]/.NDSolve[{y''[t]+((ω0)^2) y[t]-λ ((y[t])^3)==0, y[0]==A,y'[0]==0},y[t],{t,0,T},PrecisionGoal->∞][[1]];
ω0
wA
Plot[{xN[t],xA,A Cos[ω0 t]},{t,0,T},PlotStyle->{{Black,Thickness[0.01]},{Red,Thickness[0.008],Dashed},{Green,Thickness[0.004]}}]

The perturbative solution still does a good job in describing the oscillator while the harmonic approximation clearly breaks down.  If we push too hard, our perturbative solution will also break down.

A=1.8;
xN[t_]=y[t]/.NDSolve[{y''[t]+((ω0)^2) y[t]-λ ((y[t])^3)==0, y[0]==A,y'[0]==0},y[t],{t,0,T},PrecisionGoal->∞][[1]];
ω0
wA
Plot[{xN[t],xA,A Cos[ω0 t]},{t,0,T},PlotStyle->{{Black,Thickness[0.01]},{Red,Thickness[0.008],Dashed},{Green,Thickness[0.004]}}]