Week 3 - Solving second order ODEs

Aim:- To write a program that solves the following ODE.

`(d^2theta)/(dt^2)+b/m (d theta)/dt+g/Lsintheta=0`

To write a program in Matlab that will simulate the pendulum motion.

Explanation:-

In the above equation,

g = gravity in m/s2 = 9.81

L = length of the pendulum in m = 1

m = mass of the ball in kg = 1

b=damping coefficient = 0.05

To solve the higher order differential equation we need to convert it into lower order differential equation so from above equation we have;

`theta=theta1`

`(d theta)/dt=(d theta1)/dt=theta2`

`(d^2theta1)/(dt^2)=d/dt((d theta1)/dt)=(d theta2)/dt`

`(d^2theta1)/(dt^2)+b/m(d theta1)/dt+g/Lsintheta=0`

`(d theta2)/dt+b/m(theta2)+g/Lsintheta1=0``(d theta2)/dt=-b/mtheta2-g/Lsin theta1`

Matlab Code:-

Matlab Code:-

clear all
close all
clc
 

b = 0.05;
g = 9.81;
l = 1;
m =1;
 
%initial condition
theta_0 = [0;3];
 
%time points
t_span = linspace(0,20,600);
 
% solving 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));
ylabel("Amplitude");
xlabel("Time");
axis([0 20 -6 6]);
grid on;
 
%code for motion of pendulum
figure (2)
ct = 1;
for i=1:length(results(:,1))
    x0=0;
    y0=0;
    plot([-1 1],[0 0],'linewidth',5,'color','m');
    hold on
    %position of free end
    x1=l*sin((results(i,1)));
    y1=-l*cos((results(i,1)));
    
    %plotting the length of pendulum
    plot([x0,x1],[y0,y1],'linewidth',3,'color','g');
    
    %ploting of fixed end
    plot(x0,y0,'s','MarkerFaceColor','k','MarkerSize',12);
   
    %plotting of bob
    plot(x1,y1,'o','MarkerFaceColor','y','MarkerSize',25);
    hold off
    axis([-1.5 1.5 -1.6 0.5]);
    grid on
    titlestring=['Pendulum motion at time  t=’,num2str(t(i)),'seconds'];
    title(titlestring,'fontsize',8);
    
    M(ct)=getframe(gcf);
    ct=ct+1;
end


%Code for animation
movie(M)
videofile = VideoWriter('Pendulum motion.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

ODE is written below using variable t,theta,b,g,L,m  in ode_func Function.This function will be called to solve ODE .

% Code for ode_func
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

Result:-

1) A Matlab code is written to solve ODE and position of pendulum w.r.t time is found out which is attached here, Row 1 shows time and Row 2 shows Position of pendulum.

2)The simulation of Pendulum moton is done, which is given in following link:-

https://drive.google.com/file/d/1B84JbtH-bNqVdPS6r4CftrJ4yCdbe0SU/view?usp=sharing