Week 3 - Solving second order ODEs

Aim: write a code to get 2ode oscillation of a simple pendulum

Objective: we are going to sort the 2ode into 1ode by and by using mat lab coding methods and finally we are going to get a video file as a result which shows an oscillation of a simple pendulum

Theory: a simple pendulum is on end is fixed point and other end is free so if we oscillate a pendulum with some initial velocity i.e. 3m/s the displacement of a bob will changes with respect to time

basically, a simple pendulum involved with 2ode one is velocity and another one is displacement with respect to time

In a mathematical form we can write is as

d^2Ѳ/dt^2+(b/m*dѲ/dt) +g/l*sinѲ=0

where l=1; %length of the string mi m

b=0.05; %damping kg/m

g=9.81; %gravity in m/s^3

m=0.1; %mass in kg

Ѳ=theta in degrees

further we need to break it into 1ode

so, if Ѳ= Ѳ1

dѲ/dt= dѲ1/dt= Ѳ2

d^2Ѳ/dt^2= d^2Ѳ1/dt^2==(dѲ1/dt) == dѲ2/dt

dѲ2/dt+ b/m* Ѳ2+g/l*sin Ѳ1=0

so dѲ1/dt= Ѳ2; dѲ2/dt= -b/m Ѳ2-g/l*sin Ѳ1

up to now we got mathematical solution then we are going to code a program in mat lab.

Code and explanations:

as a basic thing we must clear the editor window, by using clear all, close all, clc

firstly, we need to give input values of a simple pendulum

such as l=1; %length of the string im m

b=0.05; %damping kg/m

g=9.81; %gravity in m/s^3

m=0.1; %mass in kg

Ѳ=theta in degrees

And we are going to describe the velocity and displacement values with respect to theta as given

Then we are giving the t_span values as required by using linspace command

And next the most important thing we are going to solve a 2ode by using an ode45 command, in that we need to describe the ode45 inputs such as challenge,2^nd ode, t_span, theta conditions, then only ode45 command will get execute.

In addition to that we need to call a function command for 2^nd ode in function

Finally, we will get one of our output by using plot commands for velocity and displacement with respect to time

The for-2nd plot we are going to get the position radius with respect to angular velocity of the bob

And most important thing we are going to create an oscillation of a pendulum by using a for loop command and, we are creating a bod by using a markers command.

Here by using a trigonometry we are going to determine the oscillation of a pendulum against theta so for number of theta values we will get a different position of a pendulum from this for loop command

So, all these are going to be captured with movie command.

Conclusion: Finally, we will get a video as expected.

Errors faced:

I was faced an error in for loop i.e. wrong calculation

Through support session I rectified it, by putting ‘-‘ in y1 values.

And I got confused at how subplot command works

Form support session I got the solution.

https://drive.google.com/file/d/1sGT16iFGHN1p1fRVTfW_DZXt9ne9JRiO/view?usp=sharing

% clearing the screen,command line, work space and closing all plots
clear all
close all
clc

%inputs

l=1;  %length of the string im m
b=0.05; %damping kg/m
g=9.81; %gravity in m/s^3
m=0.1;  %mass in kg

%initial conditions
% theta inputs are intial displacement and velocity

theta_0=[0,3];

%time points

t_span = linspace(0,20,100);

%[t,results]=ode45('challenge 2nd_ode',t_span,theta_0);

[t,results]=ode45(@(time,theta) pendulum_damped_ode(time,theta,b,g,l,m),t_span,theta_0);

%plotting the velocity and displacement values

subplot(5,4,[13,14,17,18]);
plot(t,results(:,1),'color','b');
hold on
plot(t,results(:,2),'color','r');
xlabel('time');
ylabel('velocity-displacement');
grid on;
hold off;

%plotting the position radius with respect to angular velocity of the bob

subplot(5,4,[15,16,19,20]);
plot(results(:,1),results(:,2));
xlabel('position_radius')
ylabel('angular-velocity')
grid on;

%writting for loop cod for a pendulum oscillation

ct=1;

for i=1:length(results(:,1))
  x0=0;
  y0=0;
  x1=l*sin(results(i,1));
  y1=-l*cos(results(i,1));
  subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);
  plot([-1 1],[0 0],'linewidth',4,'color','k');
  hold on;
  plot([x0 x1],[y0 y1],'linewidth',2,'color','b');
  hold on;
  

  % by using markers command we are goining to get a bod
  
  plot(x1,y1,'-o','markers',20,'markerfacecolor','r');  %plotting the bob using the marker and markerfacecolor
  
 
 hold on;
 axis equal;
 grid on;
 axis([-1.5 1.5 -1.5 0.5]); 
 pause(0.3);
 hold off;

 M(ct)=getframe(gcf);
 ct=ct+1;





end

 %creating a avifile by using a movie command

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


function[dtheta_dt]=pendulum_damped_ode(time,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