Week 4.1 - Solving second order ODEs

AIM : To write a code in matlab that will simulate pendulum motion

PROCEDURE :   

                      STEPS  FOLLOWED WHILE WRITING CODE 

• FIRST GIVE INPUTS OF GIVEN VALUES AND CONSTANTS
• THEN WRITE DOWN THE INITIAL CONDITION GIVEN IN CHALLENGE
• THEN GIVE THE TIME POINTS FOR MOTION AS MENTIONED UPTO 20 SEC
• THEN SOLVE ODE USING ODE45 INBUILT MATLAB FUNCTION
• FOR THAT WE NEED TO CALL THE FUNCTION WHICH WILL CALCULATE VALUES FOR ANGULAR DISPLACEMENT &ANGULAR VELOCITY
• THEN WE NEED TO DO PLOTING FOR ANGULAR VELOCITY AND ANGULAR DISPLACEMENT 
• THEN WE NEED TO ANIMATION PART
• FOR ANIMATION WE NEED TO DEFINE POINTS
• FIRST POINT IS FIXED 
• FOR SECOND POINT WE NEED DEFINE IT ASIT CHANGES WITH RESPECT TO TIME AND ANGLE
• GIVING COUNTER COMMAND COMMAND AS IT IS REPRESENTING EVERY FRAMES
• FOR SOLVING IS REPEATITVELY WE NEED TO USE FOR LOOP
• BY USING PLOT COMMAND AND GIVEN LOCATIONS  FOR SUPPORT ROD , FIXED POINT , LENGTH OF ROD
• SIMILARLY FOR BOD PENDULUM WE HAVE TO PLOT BY USING MARKER COMMAND FOR SPHERE SHAPE

 USING MOVIE COMMAND FOR ANIMATION CREATION

 CODE FOR MATLAB :

clear all
close all
clc

%Given values or consider this values as constant
b=0.05;
g=9.81;
l=1;
m=1;

%initial condition
theta_0 = [0 3];

%time points
t_span = linspace(0,20,500);

% solve ODE using ode45(matlab inbuilt function for \)
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

%plotting for tracing of angular velocity and angular displacement
plot(t,results(:,1));%angular displacement
hold on
plot(t,results(:,2));%angular velocity
ylabel('Plot')
xlabel('time')

% ANIMATION

X0=0;
Y0=0;

ct=1;% counter actually represents each and every frame 
for i=1:length(t)
    X1=l*sin(results(i,1));
    Y1=-l*cos(results(i,1));
    
    figure(2)
    plot([-1 1],[0 0],'linewidth',6,'color','r');
    axis([-2 2 -2 2]);
    line([X0 X1],[Y0 Y1],'linewidth',5,'color','y')
    hold on;
    plot(X0,Y0,'O','markersize',3,'markerfacecolor','b');
    plot(X1,Y1,'O','markersize',30,'markerfacecolor','g');
    title('pendulum motion');
    legend('support rod','length of pendulum','fixed point','bod of pendulum');
    grid on;
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
end

movie(M);
videofile=videowriter('pendulum.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M);
close(videofile);
    

FUNCTION CODE : 

%FUNCTION CODE
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)
       
    theta1 = theta(1)
    theta2 = theta(2)
    dtheta1_dt = theta2;%first element is angular velocity , which is at theta2 position
    dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);%this wil obtain after conversion of 2nd order ode to 1st order ode
    dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

PLOT FOR ANGULAR DISPLACEMENT AND VELOCITY :

Pendulum motion figure :