function pencompare(theta0,omega0,N,t1,beta)
% PENCOMPARE(theta0,omega0,N,t1,beta)
% Use a 4th order Runge-Kutta method to compare
% the numerical solutions of the simple pendulum
% and of its linearization (harmonic oscillator).
% The optional fifth argument provides for changing
% the damping from the default zero to any value.
%
% If no initial conditions theta0 and omega0 are
% specified, the default values Pi/2 and 0 are used.
%
% Comparison meaningful for small angles.  Not really 
% a surprise: Try out pencompare(0,2.1,1000,50,0.1)
% and pencompare(0,2.05,1000,50,0.1).
% The solution is correct! Explain what is going on.
%
% The values of the optional parameters N and t1
% determine the number of steps taken, and the 
% final time. The default values are 1000 and 20. 
%
% All rights reserved. Matt Kawski. April 1999.
% http://math.la.asu.edu/~kawski
%
if nargin < 5, beta=0; end
if nargin < 4, t1=20; end
if nargin < 3, N=1000; end
if nargin < 2, omega0=0; end
if nargin < 1, theta0=pi/12; end
%
% Define nominal values of the system parameters
% mass of the pendulum, length, grav accel, damping
m=1
L=9.81
g=9.81
%
% Define the time interval on which to numerically 
% solve the DE. Specify the number of subintervals 
% to be used. 
% For efficiency reserve contiguous memory blocks 
% to store the values generated by the integration 
% algorithm. Define the first value in each array. 
%
t0=0
dt=(t1-t0)/N
tt=[t0:dt:t1];
ttheta=zeros(size(tt));
oomega=zeros(size(tt));
ttheta(1)=theta0;
oomega(1)=omega0;
ltheta=ttheta;
lomega=oomega;
%
% Use a 4th order Runge Kutta method to solve the system of DEs.
%
for k=1:N
   a1=dt*(-beta/m/L*oomega(k)-g/L*sin(ttheta(k)));
   b1=dt*oomega(k);
   a2=dt*(-beta/m/L*(oomega(k)+a1/2)-g/L*(sin(ttheta(k)+b1/2)));
   b2=dt*(oomega(k)+a1/2);
   a3=dt*(-beta/m/L*(oomega(k)+a2/2)-g/L*(sin(ttheta(k)+b2/2)));
   b3=dt*(oomega(k)+a2/2);
   a4=dt*(-beta/m/L*(oomega(k)+a3/2)-g/L*(sin(ttheta(k)+b3/2)));
   b4=dt*(oomega(k)+a3);
   oomega(k+1)=oomega(k)+1/6*(a1+2*a2+2*a3+a4);
   ttheta(k+1)=ttheta(k)+1/6*(b1+2*b2+2*b3+b4);
   a1=dt*(-beta/m/L*lomega(k)-g/L*ltheta(k));
   b1=dt*lomega(k);
   a2=dt*(-beta/m/L*(lomega(k)+a1/2)-g/L*(ltheta(k)+b1/2));
   b2=dt*(lomega(k)+a1/2);
   a3=dt*(-beta/m/L*(lomega(k)+a2/2)-g/L*(ltheta(k)+b2/2));
   b3=dt*(lomega(k)+a2/2);
   a4=dt*(-beta/m/L*(lomega(k)+a3/2)-g/L*(ltheta(k)+b3/2));
   b4=dt*(lomega(k)+a3);
   lomega(k+1)=lomega(k)+1/6*(a1+2*a2+2*a3+a4);
   ltheta(k+1)=ltheta(k)+1/6*(b1+2*b2+2*b3+b4);
end
%
% Plot the solutions.
%
clf
hold on
plot(tt,oomega,'b-',tt,ttheta,'r-')
plot(tt,lomega,'c-',tt,ltheta,'m-')
title('Numerical comparison of simple and linearized pendulum')
ww=axis;
w4=ww(4)
text(t1/10,0.9*w4,'Blue/red simple. Cyan/magenta=linearized')

%
% compare with the the exact solutions of the linearized pendulum
% 
p=beta/m/L;
q=g/L;
if 4*q-p^2 <= 0 then
text(t1/3,0.7*max(ttheta),'Critically or overdamped pendulum')
else
   omega=sqrt(q-p^2/4)
   A=theta0
   B=omega0+theta0*p/2
exactth=(A*cos(omega*tt)+B*sin(omega*tt)).*exp(-p/2*tt);
exactom=((B-A*p/2)*cos(omega*tt)+...
         (-A-B*p/2)*sin(omega*tt)).*exp(-p/2*tt); 
plot(tt,exactom,'g-',tt,exactth,'y-')
plot(tt,lomega,'c-',tt,ltheta,'m-')
text(t1/10,0.8*w4,'Exact soln of linearized: green/yellow')
end