(**************************************)
(*Initialisation*)(* 1D SHM with Euler method *)
Deltat=0.01;tini=0;l=1.0;g=9.81;
omega[0]=0;theta[0]=2*(Pi/180);
T=(2Pi)Sqrt[l/g];
tlast=20T;
(**************************************)

(**************************************)
(* Calculation *)
ksq=g/l;
Ni=IntegerPart[(tlast-tini)/Deltat];
For[i=0,i≤Ni,i++,
t[i]=tini+i*Deltat;
omega[i+1]=omega[i]+Deltat*(-ksq*theta[i]);

theta[i+1]=theta[i]+Deltat*omega[i];

];(*end for *)
tttheta= Table[{t[i],theta[i]},{i,0,Ni}];
ttomega= Table[{t[i],omega[i]},{i,0,Ni}];



(**************************************)

(* Output *)
(**************************************)
ListPlot[tttheta,PlotLabel->"theta vs time "]
ListPlot[ttomega,PlotLabel->"omega vs time "]

(**************************************)

(**************************************)
(*Initialisation*)(* 1D SHM with Euler-Cromer method *)
Deltat=0.01;tini=0;l=1.0;g=9.81;
omega[0]=0;theta[0]=2*(Pi/180);
T=(2Pi)Sqrt[l/g];
tlast=20T;
(**************************************)

(**************************************)
(* Calculation *)
ksq=g/l;
Ni=IntegerPart[(tlast-tini)/Deltat];
For[i=0,i≤Ni,i++,
t[i]=tini+i*Deltat;
omega[i+1]=omega[i]+Deltat*(-ksq*theta[i]);

theta[i+1]=theta[i]+Deltat*omega[i+1];

];(*end for *)
tttheta= Table[{t[i],theta[i]},{i,0,Ni}];
ttomega= Table[{t[i],omega[i]},{i,0,Ni}];



(**************************************)

(* Output *)
(**************************************)
ListPlot[tttheta,PlotLabel->"theta vs time "]
ListPlot[ttomega,PlotLabel->"omega vs time "]

(**************************************)

(**************************************)
(*Initialisation*)(* 1D SHM, with Euler-Cromer method, with damping *)
Deltat=0.01;tini=0;l=1.0;g=9.81;
omega[0]=0;theta[0]=2*(Pi/180);
T=(2Pi)Sqrt[l/g];m=1.0;k=0.350;
q=k/m;
tlast=5T;
(**************************************)

(**************************************)
(* Calculation *)
ksq=g/l;
Ni=IntegerPart[(tlast-tini)/Deltat];
For[i=0,i≤Ni,i++,
t[i]=tini+i*Deltat;
omega[i+1]=omega[i]+Deltat*(-ksq*theta[i])
-Deltat*(omega[i]*q)
;
theta[i+1]=theta[i]+Deltat*omega[i+1];

];(*end for *)
tttheta= Table[{t[i],theta[i]},{i,0,Ni}];
ttomega= Table[{t[i],omega[i]},{i,0,Ni}];



(**************************************)

(* Output *)
(**************************************)
ListPlot[tttheta,PlotLabel->"theta vs time ",PlotRange->{-theta[0],theta[0]}]
ListPlot[ttomega,PlotLabel->"omega vs time ",PlotRange->{-theta[0]*Sqrt[g/l],theta[0]*Sqrt[g/l]}]

(**************************************)

(**************************************)
(*Initialisation*)(* 1D SHM, with Euler-Cromer method, with damping and driving force *)
Deltat=0.01;tini=0;l=1.0;g=9.81;
omega[0]=0;theta[0]=2*(Pi/180);
T=(2Pi)Sqrt[l/g];m=1.0;k=0.0;
q=k/m;FD=1.0;OmegaD=0.50Sqrt[g/l];
tlast=100T;
(**************************************)

(**************************************)
(* Calculation *)
ksq=g/l;
Ni=IntegerPart[(tlast-tini)/Deltat];
largesttheta=0;

For[i=0,i≤Ni,i++,
t[i]=tini+i*Deltat;
omega[i+1]=omega[i]+Deltat*(-ksq*theta[i])
-Deltat*(omega[i]*q)+Deltat*(FD/m)*Sin[OmegaD*i*Deltat];
theta[i+1]=theta[i]+Deltat*omega[i+1];

largesttheta= Max[largesttheta,Abs[theta[i+1]]];

];(*end for *)
tttheta= Table[{t[i],theta[i]},{i,0,Ni}];
ttomega= Table[{t[i],omega[i]},{i,0,Ni}];

Print["largesttheta= ",largesttheta];

(**************************************)

(* Output *)
(**************************************)
Amp=2.0(FD/m)/Sqrt[((g/l)-OmegaD^2)^2+(q*OmegaD)^2];

ListPlot[tttheta,PlotLabel->"theta vs time ",PlotRange->{-largesttheta,largesttheta}]
ListPlot[ttomega,PlotLabel->"omega vs time ",PlotRange->{-largesttheta*Sqrt[g/l],largesttheta*Sqrt[g/l]}]

(**************************************)