Week 3 - Solving second order ODEs

1. Program for simulation of simple pendulum :

clear all
close all
clc

% input 
b = 0.05; % damping coefficient
l = 1;    % pendulam rod length
g = 9.81; % acceleration due to gravity
m = 1;    % pendulam bob mass

% initial conditionb
theta_0=[0;3];

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


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

%plots for tracing the angular displacment and angular velocity
figure(1)
plot(t,results(:,1),'LineWidth',3,'Color','b');
hold on
plot(t,results(:,2),'LineWidth',3,'Color','r');
xlabel('Time t in sec');
ylabel('Plot');
grid on

% Animation of Single bob preforming simple pendulam motion

X0 = 0;
Y0 = 0;

ct = 1
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',5,'Color','black');
    axis([-2 2 -2 2]);
    line([X0 X1],[Y0 Y1],'linewidth',2,'color','g');
    hold on;
    plot(X0,Y0,'O','MarkerSize',4,'markerfacecolor','b');
    plot(X1,Y1,'O','MarkerSize',20,'markerfaceColor','r');
    legend('support','length of pendulum','fixed point','bob of pendulum');
    grid on;
    hold off;
    M(ct)=getframe(gcf)
    ct = ct + 1
end

movie(M);
videofile=VideoWriter('motion_pendulum.avi','uncompressed_file avi');
open('videofile');
VideoWrite(videofile,M);
close('videofile')

2. Program code of function which we are going to call :

function [dtheta_dt] = ode_func(t,theta,b,l,g,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

Explaination :

Stepwise explaination of program as follows

1. In the initial , i have written the given values of damping coefficient , length of pendulum rod, acceleration due to gravity and mass of the pendulum bob as you can see in the code .

2.Then initial condition written by the command : theta_0 = [0;3] : which means the value of theta is givan for initial displacement is 0 and the initial velocity is 3 at time t= 0.

3.Using the command : t_span = linspace(0,20,400) : which basically means the 400 time points of movement of bob pendulum created , bob pendulum starting from 0 and ending after 20 sec.

4.Then ODE solver (Command : ode(@(t,theta) diff_func(t,theta,b,l,g,m),t_span,theta_0)) function is used to solve the 2nd order derivatives which is the function which shows the movement of bob pendulum

5. Code for ode function is above's 2nd code , after exicuting the function there was an error can be seen below attached image name 1 ode , error is resolved just changing the name where error is found.

6.For tracing the Angular displacement and the Angular Velocity i have used the command :plot(t,results(:,1),'linewidth',3,'color','b') : which plot the graph of Angular Displacement with respect to time as this

command contain the term (t,results(:,1)) : because the 1st column of results is column of Angular Displacement and same we can plot the graph of Angular Velocity just by changing the term (t,results(:,2)).

7.Labels are provided and after exicuting the code we got the plot can be seen the below attached image name

8. After exicuting the code the results we got can be seen in below attached image name 2 ode

9. Then to make it into the animation the coordinates difined as X0,X1,Y0,Y1

10. For loop command is used :for i = 1:length(t) : which means the loop(action) will start from 1 and will goes on till the time t defined and creat the position of bob pendulum for every value of time and displacement .

11. new coordinates created by for loop written as X1=l*sin(results(i,1))and Y1=-l*cos(results(i,1))

12. Then Plot is plotted using the plot command : plot([-1 1],[0 0],'linewidth',5,'Color','black') Command for Support.
axis([-2 2 -2 2]) : Command for Axis limites
line([X0 X1],[Y0 Y1],'linewidth',2,'color','g') : For showing the bob pendulum rod
hold on; : To hold the every plot created
plot(X0,Y0,'O','MarkerSize',4,'markerfacecolor','b'): For locating the fixed point
plot(X1,Y1,'O','MarkerSize',20,'markerfaceColor','r'): For locateing the bob with size
legend('support','length of pendulum','fixed point','bob of pendulum') : this command is used for showing the name of the componants

13. For making the movie of bob pendulum the counter command is used as : ct = ct+1 : which creates the counter the every frame .

14. Movie command [M(ct)=getframe(gcf)] is used to create the movie of the every frame for every position ,we have created of pendulum.

15. Exicuting the program there was the error can be seen in below attached image name 4 and encounterd [ movie(ct) ] this error just rewritting the command movie(M) .

16.After ending the program we needed the file where we need to dumb this file created the file outside the for loop as :

movie(M)
videofile=VideoWriter('motion_pendulum.avi','uncompressed_file avi')
open('videofile')
VideoWrite(videofile,M)
close('videofile')

14. After Executing the code we got the animation of pendulum performing the motion which can be seen below attached file name (2021-12…8-00.mp4) or by clicking the given link : https://youtu.be/JkufNf_12KA