function pendulum.m
function uprime = pendulum(t,u)
%Pendulum(t,x) returns the state derivatives
%of the simple pendulum equation
theta = u(1); % u is the state vector (theta and derivs)
thetadot = u(2);
g = 9.8; % acceleration due to gravity
L = 1; % length of pendulum is 1 meter
uprime(1) = thetadot; % define u1' = u2 = thetadot
uprime(2) = -(g/L) * sin(theta); % define u2'
script file swing.m
initial=[5/180*pi,0]';
%initial conditions 0 velocity at 5 degrees
[t,numtheta] = ode23('pendulum',0,20,initial);
subplot(211),plot(t,numtheta(:,1)),...
title('theta'),xlabel('time(s)'),...
ylabel('radians'),grid,...
subplot(212),plot(t,numtheta(:,2)),...
title('theta dot'),xlabel('time(s)'),...
ylabel('radians'),grid
Challenge solution (put in file SwingPlus.m):
clear
x0 = input('Enter initial angle in x-direction (degrees) ');
x0 = x0/180*pi;
y0 = input('Enter intial angle in y-direction (degrees) ');
y0 = y0/180*pi;
vx0 = input('Enter intial velocity in x-direction (m/s) ');
vy0 = input('Enter initial velocity in y-direction (m/s) ');
initial_x=[x0,vx0]'; %initial conditions in x
initial_y=[y0,vy0]'; %initial conditions in y
[time_x,numtheta] = ode23('pendulum',0,20,initial_x);
[time_y,numphi] = ode23('pendulum',0,20,initial_y);
subplot(2,2,1),plot(time_x,numtheta(:,1)),...
title('theta'),xlabel('time(s)'),ylabel('radians'),grid,...
subplot(2,2,3),plot(time_x,numtheta(:,2)),...
title('theta dot'),xlabel('time(s)'),...
ylabel('radians/s'),grid,...
subplot(2,2,2),plot(time_y,numphi(:,1)),...
title('phi'),xlabel('time(s)'),...
ylabel('radians'),grid,...
subplot(2,2,4),plot(time_y,numphi(:,2)),...
title('phi dot'),xlabel('time(s)'),...
ylabel('radians/s'),grid
pause;
nx = length(time_x);
ny = length(time_y);
if (nx > ny)
shorttx = time_x(1:ny);
shortty = time_y;
shortnumtheta = numtheta(1:ny,1);
shortnumphi = numphi(:,1);
elseif (ny > nx)
shorttx = time_x;
shortty = time_y(1:nx);
shortnumtheta = numtheta(:,1);
shortnumphi = numphi(1:nx,1);
else
shorttx = time_x;
shortty = time_y;
shortnumtheta = numtheta(:,1);
shortnumphi = numphi(:,1);
end
clg
plot(cos(shortnumtheta),sin(shortnumphi))
%axis([-1 1 -1 1])