Week 3 - Solving second order ODEs

clear all
close all
clc


b=0.05; % damping coefficient
g=9.81; %in m/s^2
l=1; % in meter
m=1; % in kg


% initial condition
theta_0=[0;3];

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

% solve ODE
[t,results]= ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))

xlabel('time in sec')
ylabel('angular velocity(radian/sec) and displacement(radian)')

clear all
close all
clc

% inputs
b=0.05;
g=9.81; %in m/s^2
l=1; % in meter
m=1; % in kg

ct=1;
% initial condition
theta_0=[0;3];

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

% solve ODE
[t,results]= ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
hold off
xlabel('time in sec')
ylabel('angular velocity(radian/sec)and displacement(radian)')

for i=1:length(t_span)
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));

    %supporting bar
plot([-1 1],[0 0],'linewidth',2,'color','r')
hold on
% pendulum length
plot([0 x1],[0 y1],'linewidth',3,'color','b')
% Bob 
plot(x1,y1,'o','MarkerSize',20,'MarkerFaceColor','r')
axis([-4 4 -3 1])

pause(0.3);
M(ct)=getframe(gcf);
hold off
ct=ct+1;
end

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

link for video

after download the video you be enable to see the video

https://1drv.ms/v/s!AjWNeQaHir-bgUtWb-LtRcjJb296?e=P2fYto

function  [dtheta_dt]= ode_func(t,theta,b,g,l,m)


theta1= theta(1);
theta2= theta(2);
dtheta1_dt= theta2;
dtheta2_dt= -(b/m)*theta2 -(g/l)*sin(theta1);
dtheta_dt= [dtheta1_dt; dtheta2_dt];
end

error1

In the function i am not using ')' bracket after the sin(theta1; 

function  [dtheta_dt]= ode_func(t,theta,b,g,l,m)


theta1= theta(1);
theta2= theta(2);
dtheta1_dt= theta2;
dtheta2_dt= -(b/m)*theta2 -(g/l)*sin(theta1;
dtheta_dt= [dtheta1_dt; dtheta2_dt];
end

error2

In pendulum length i will use the co-ordinate 
plot([0 x1],[0 -y1],'linewidth',3,'color','b')

the rotation will be opposite as desire

error3

first i have to use 

x1=l*sind(i);
y1=l*cosd(i);

the bob revolve above the support

error4

then i have to co-ordinate as

x1=l*sind(i,1);
y1=l*cosd(i,1);

it will show error using sind

error5

I am not able to make bob
so I have to use * as bob

then i know abot markersize which clear show the bob and its size