%Solve the large angle pendulum problem using rk4 % Izzi Urieli - December 8, 2007 clear all THETA = 1; % index of angle [radians] OMEGA = 2; % index of angular velocity [rads/sec] HEIGHT = 3; % index of pendulum mass height [m] global g length g = 9.807; % [m/s/s} accel. due to gravity length = input('Enter pendulum length [m]: '); angle0 = input('Enter initial pendulum angle [degrees]: '); period = 2*pi*sqrt(length/g); fprintf('small angle period is %.3f seconds\n',period); tfinal = 1.5*period; %######################using rk4: n = input('Enter the number of solution points: '); dt = tfinal/(n - 1); t = 0; y(THETA) = angle0*pi/180; % initial pendulum angle in radians y(OMEGA) = 0; y(HEIGHT) = length*(1 - cos(y(THETA))); % initial pend. height time(1) = t; angle(1) = angle0; height(1) = y(HEIGHT)*100; % cm (for scaling) for i = 2:n [t, y, dy] = rk4('dpend',2,t,dt,y); time(i) = t; angle(i) = y(THETA)*180/pi; height(i) = y(HEIGHT)*100; end %##############################plotting: plot(time,angle) grid on hold on axis([0 tfinal -90 90]) x0 = [0,tfinal]; y0 = [0,0]; plot(x0,y0,'k-') x1 = [period,period]; y1 = [-90,90]; plot(x1,y1,'r-') plot(time,height, 'g-') xlabel('time [seconds]') ylabel('pendulum angle[degrees], height[cm] (green)') title('pendulum angle and height vs time')