Solution for Nonlinear Pendulum lab

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)

g = 9.8; % acceleration due to gravity

L = 1; % length of pendulum is 1 meter

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)'),...

subplot(212),plot(t,numtheta(:,2)),...

title('theta dot'),xlabel('time(s)'),...

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)),...

subplot(2,2,3),plot(time_x,numtheta(:,2)),...

title('theta dot'),xlabel('time(s)'),...

subplot(2,2,2),plot(time_y,numphi(:,1)),...

title('phi'),xlabel('time(s)'),...

subplot(2,2,4),plot(time_y,numphi(:,2)),...

title('phi dot'),xlabel('time(s)'),...

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])