Week 3 - Solving second order ODEs

Here Solving the ODE in the sense , when you have ODE , break it down to first order ODE , we have second Order ODE in this challenge.

Solving theta ( Displacement ) and also solving rate at which the theta is changing with respect to time , which is nothing but ( velocity ) .

Here the given ODE represents the equation of motion of a simple pendulum with damping

In the above equation ,

g= Gravity in m/s^2

L=Length of simple pendulum in m.

b=Damping co-efficient.

m=mass of the ball in kg

Given values are g=9.81m/s^2 ; L=1 m ; b=0.05 ; m= 1 kg.

the following code will give the position of pendulum at various position with respective to time and angles( dispalcement)

clear all
close all
clc

b = 0.05;
g = 9.81;
l = 1;
m = 1 ;

%initial conditions
theta_0 = [0;5];

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

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

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

ylabel('plot')
xlabel('time')

fig=1 ;
for i=1:length(t)
    
   %intial position of Pendulum
   
   x0=0;
   y0=0;
   %final position of Pendulum
   %here we will find the position of pendulum at various points w.r.t time and angles
   
   x1= 1*sin(results(i,1));
   y1= -1*cos(results(i,1));
   
   %Plotting as
   figure(2)
   
   %baseline
   plot([-1,1],[0,0],'linewidth',6,'color','k');
   
   %Pendulum line
   line([x0,x1],[y0,y1],'linewidth',3,'color','b');
   hold on
   %pendulum ball
   plot(x1,y1,'o','markers',25,'markerfacecolor','g');
   grid on 
   axis([-2 2 -2 2]);
   pause(0.01);
   hold off;
   
   M(fig) = getframe(gcf);
   fig = fig+1;
end

% Storing the animation video file in avi format
movie(M)
videofile=VideoWriter('second_order_equation_for_pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Here the function code is as the following , where we will be simplifying Second order ODE into two first order ODE eqautions.

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

OUTPUT :

Variables output :

Youtube Link for animation video for position of pendulum : https://youtu.be/JypjtHLq0ws