a new title
Wrote a Matlab script to solve the non-linear pendulum equation and simulate it. The equation being solved is:
where b is the damping coefficient (0.05), m is the mass of the pendulum bob (1 kg), L is the length of the pendulum (1 m), theta is the angular displacement and g is the gravitational acceleration (9.82 m/s2). The equation is solved numerically using the ODE45 function with an initial angular velocity of 3 rad/s and angular displacement of 0 radians for 20 sec.
To view the animation please visit:
https://www.youtube.com/watch?v=CByQidAC1V4&feature=youtu.be
clear allclose allclcg = 9.81; %Gravitational accelaration (ms^-2)L = 1; %Length of the pendulum (m)m = 1; %Mass of the pendulum (kg)b = 0.05; %Damping coefficient runtime = [0 20]; %Simulating b/w 0 and 20sinicondi =[0; 3]; %Initial condition @ 0 angular displacement and ...angular veliocity of 3 rad/sectheta = @(t, y) [y(2); (-(g/L)*sin(y(1))-(b/m)*y(2))]; %Fucntion handle for Pendulum Equation[t, sol] = ode45(theta, runtime, inicondi); x_cord = L*sin(sol(:,1)); %Cartesian Coordinates from angular-displacememnty_cord = -L*cos(sol(:,1));x_vel = L.*sol(:,2).*cos(sol(:,1));y_vel = L.*sol(:,2).*sin(sol(:,1));figure(1)%===Plot wall===%%Horizontal linesplot([-1.5 1.5],[0 0],'k','linewidth',2);hold onplot([-1.5 1.5],[0.125 0.125],'k','linewidth',1);plot([-1.5 1.5],[0.25 0.25],'k','linewidth',1);plot([-1.5 1.5],[0.375 0.375],'k','linewidth',1);plot([-1.5 1.5],[0.5 0.5],'k','linewidth',2);%Vertical linesplot([-1.2 -1.2],[0 0.125],'k','linewidth',1);for lev = 0:3 for brk = 0:10 tot = 0; if(lev==0 || lev==2) if(tot > -1.5 && tot <= 1.5) plot([-1.5 -1.5]+(0.3*brk),[0 0.125]+(0.125*lev),'k','linewidth',1); tot =+ 0.3; end else if(tot > -1.5 && tot < 1.5) plot([-1.35 -1.35]+(0.3*brk),[0 0.125]+(0.125*lev),'k','linewidth',1); tot =+ 0.3; end end endend %===Plot Pendulum Swing Axis===%plot (0, 0, 'sk','markersize',10);%===Plot Pendulum Vertical Reference Axis===%plot([0 0], [0.75 -1.25],'-.k');xlim([-1.5 1.5])ylim([-1.5 1.5])for i = 1:length(t) figure(1);%===Plot 1m rod===%p1 = plot([0 x_cord(i)],[0 y_cord(i)],'b','LineWidth',3);%axis off;%===Plot Pendulum bob===%p2 = plot(x_cord(i),y_cord(i),'o','markersize',16,'markerfacecolor','c','MarkerEdgeColor','k');%===Velocity vector===%p3 = quiver(x_cord(i), y_cord(i), 0.1*x_vel(i), 0.1*y_vel(i), '--r','maxheadsize',1);%legend(p3,'Linear Velocity Vector');%===Pendulum Bob Tracer===%h1 = plot(x_cord(i), y_cord(i),'.k','markersize', 0.2);drawnow;pause(0.002); %%=====Angular Displacement & Angular Velocity Plots=====%%%===Ploting using subplot===%% figure(2);% subplot(2,1,1)% plot(t(i),sol(i,1),'ok');% hold on% xlabel('Time')% ylabel('Angular Displacement')% subplot(2,1,2);% plot(t(i),sol(i,2),'or');% xlabel('Time')% ylabel('Angular Velocity')% drawnow;% hold on%===Ploting using plot===%%figure(2);%plot(t(i), sol(i,2),'ok',t(i), sol(i,1),'or')% drawnow;% hold on%===Ploting using plotyy===%%figure(2);% [hAx, hline1, hline2] = plotyy(t(i), sol(i,1),t(i), sol(i,2));% set(hAx(1),'xlim',[round(min(t(:))),round(max(t(:)))]);% set(hAx(2),'xlim',[round(min(t(:))),round(max(t(:)))]);% set(hAx(1),'ylim',[round(min(sol(:,1))),round(max(sol(:,1)))]);% set(hAx(2),'ylim',[round(min(sol(:,2))),round(max(sol(:,2)))]);% hline1.Marker = '*';% hline1.MarkerEdgeColor = 'b';% hline2.Marker = 'o';% hline2.MarkerEdgeColor = 'r';% drawnow;% hold onF(i) = getframe(gcf);if i<=(length(t)-1)delete([p1 p2 p3]); endend%%===Save the animation in avi format===%%% video = VideoWriter('Pendulum.avi','Uncompressed AVI');% video.FrameRate = 15;% open(video);% writeVideo(video,F);% close(video);