Chapter 3 
Oscillatory Motion

Sample code 3.1.1

Deltat=0.04 (*in unit of second *);
theta0=10*(Pi/180) (* in unit of radian *);
l=1 (* in unit of m *);
omega0=0.0 (*in unit of rad/second *);
NStep =IntegerPart[ (10/Deltat)+1]; (* Number of total time step *)
g=9.81;(* in unit of m/s^2 *)

omega[0]=omega0;t[0]=0;
theta[0]=theta0;

For [i=0,i≤NStep,i++,
omega[i+1]=omega[i] -(g/l)*theta[i]*Deltat;
theta[i+1]=theta[i] + omega[i]*Deltat;
t[i+1]=t[i]+Deltat;
];(* end for *)
ttheta=Table[{t[i],theta[i]},{i,0,NStep}];
ListPlot[ttheta,PlotRange->{-2,2},Joined->True,AxesLabel->{"time(s)","θ(radian)"}]

Sample code 3.2.2:

Deltat=0.04 (*in unit of second *);
theta0=10*(Pi/180) (* in unit of radian *);
l=1 (* in unit of m *);
omega0=0.0 (*in unit of rad/second *);
NStep =IntegerPart[ (10/Deltat)+1]; (* Number of total time step *)
g=9.81;(* in unit of m/s^2 *)

omega[0]=omega0;t[0]=0;
theta[0]=theta0;

For [i=0,i≤NStep,i++,
omega[i+1]=omega[i] -(g/l)*theta[i]*Deltat;
theta[i+1]=theta[i] + omega[i+1]*Deltat;
t[i+1]=t[i]+Deltat;
];(* end for *)
ttheta=Table[{t[i],theta[i]},{i,0,NStep}];
ListPlot[ttheta,PlotRange->{-0.3,0.3},Joined->True,AxesLabel->{"time(s)","θ(radian)"}]

Sample code 3.2.3: Euler-Cromer method for forced oscillator

Deltat=0.04 (*in unit of second *);
theta0=11*(Pi/180) (* in unit of radian *);
l=1 (* in unit of m *);
omega0=0.0 (*in unit of rad/second *);
NStep =IntegerPart[ (10/Deltat)+1]; (* Number of total time step *)
g=9.81;(* in unit of m/s^2 *)
omega[0]=omega0;t[0]=0;
theta[0]=theta0;

q=1(*dimensionless*);
For [i=0,i≤NStep,i++,
omega[i+1]=omega[i] -((g/l)*theta[i]+q*omega[i])*Deltat;
theta[i+1]=theta[i] + omega[i+1]*Deltat;
t[i+1]=t[i]+Deltat;
];(* end for *)
ttheta1=Table[{t[i],theta[i]},{i,0,NStep}];


q=5(*dimensionless*);
For [i=0,i≤NStep,i++,
omega[i+1]=omega[i] -((g/l)*theta[i]+q*omega[i])*Deltat;
theta[i+1]=theta[i] + omega[i+1]*Deltat;
t[i+1]=t[i]+Deltat;
];(* end for *)
ttheta5=Table[{t[i],theta[i]},{i,0,NStep}];


q=10(*dimensionless*);
For [i=0,i≤NStep,i++,
omega[i+1]=omega[i] -((g/l)*theta[i]+q*omega[i])*Deltat;
theta[i+1]=theta[i] + omega[i+1]*Deltat;
t[i+1]=t[i]+Deltat;
];(* end for *)
ttheta10=Table[{t[i],theta[i]},{i,0,NStep}];


ListPlot[{ttheta1,ttheta5,ttheta10},PlotRange->{-0.2,0.2},Joined->True,AxesLabel->{"time(s)","θ(radian)"}]

Sample code 3.2.4: Forced oscllator solved numerically with NDSolve.

g=9.81;l=1;
s[q_]:=NDSolve[{theta''[t]==-(g/l)theta[t]-q *theta'[t],theta[0]==0.2,theta'[0]==0},theta,{t,0,30}];
(*Plot[Evaluate[theta[t]/.s[1]],{t,0,10},PlotRange->{-0.2,0.2}]
Plot[Evaluate[theta[t]/.s[5]],{t,0,10},PlotRange->{-0.2,0.2}]
Plot[Evaluate[theta[t]/.s[10]],{t,0,10},PlotRange->{-0.2,0.2}]*)
Plot[Evaluate[{theta[t]/.s[1],theta[t]/.s[5], theta[t]/.s[10]}],{t,0,10},PlotRange->{-0.2,0.2}]

Sample code 3.2.5: Euler-Cromer method for forced oscillator with a driving force

Deltat=0.04 (*in unit of second *);
theta0=11*(Pi/180) (* in unit of radian *);
l=1 (* in unit of m *);
omega0=0.0 (*in unit of rad/second *);
NStep =IntegerPart[ (100/Deltat)+1]; (* Number of total time step *)
g=9.81;(* in unit of m/s^2 *)
omega[0]=omega0;t[0]=0;
theta[0]=theta0;
Omega=Sqrt[g/l];

q=1(*dimensionless*);
FDoverm=0.2;
OmegaD=0.1*Omega

For [i=0,i≤NStep,i++,
omega[i+1]=omega[i] -((g/l)*theta[i]+q*omega[i]-(FDoverm)*Sin[OmegaD*i*Deltat])*Deltat;
theta[i+1]=theta[i] + omega[i+1]*Deltat;
t[i+1]=t[i]+Deltat;
];(* end for *)
ttheta1=Table[{t[i],theta[i]},{i,0,NStep}];

ListPlot[{ttheta1},PlotRange->{-0.2,0.2},Joined->True,AxesLabel->{"time(s)","θ(radian)"}]