% Euler integration example - simple pendulum
% This is an example m-file that illustrates a direct
% implementation of an Euler integration for a set of first order ODEs.
%
% Note: The Euler method requires a small time step, and is generally
% only recommended when the states of the system undergo smooth trajectories
% with time.
%
% The potential, kinetic, and total energies are computed and plotted to 
% illustrate how energy can be used to measure errors in a simulation.
%
% 9/1/97 rgl
% 6/13/1999 revised rgl
clear; help euler_pend;
g = 9.81;
% The following computes a length that will
% make T for small oscillations equal to 1 sec
% (this is arbitrary, but helps us choose a time step).
l = g/(4*pi*pi);
T = 2*pi*sqrt(l/g)	% period of oscillation (good for small angles)
N = 1000;				% N = points per cycle
dt = T/N;				% choose time step based on T and N
n = 2;					% n = cycles
% The following defines the two state vector arrays.
x1 = zeros(1,n*N);	% x1 = angle
x2 = zeros(1,n*N);	% x2 = angular velocity
t = zeros(1,n*N);		% t = time vector
%
x1(1) = 60*pi/180;	% initial angle
x2(1) = 0;				% initial velocity (not moving initially)
t(1) = 0.0;
%
% This is the computing loop that uses the basic
% Euler algorithm to compute the next state of the system.
% 
for k = 1:n*N-1,
   t(k+1) = t(k) + dt;
   % Note, the state equations are:
   % f1 = x2
   % f2 = -g*sin(x1)/l
   f1 = x2(k) ; f2 = -g*sin(x1(k))/l;
	x1(k+1) = x1(k) + dt*f1;
	x2(k+1) = x2(k) + dt*f2;
end
%
% Compute the total energy per unit mass
%  
for k = 1:n*N,
   KE(k) = 0.5*l*l*x2(k)*x2(k);
   PE(k) = g*l*(1-cos(x1(k)));
	E(k) = g*l*(1-cos(x1(k))) + 0.5*l*l*x2(k)*x2(k);
end
%
subplot(2,1,1), plot(t,x1)
subplot(2,1,2), plot(t,PE,t,KE,t,E)
