Week 3 - Solving second order ODEs

Program to simulate the motion of a pendulum.

First of all, the required inputs are given for solving the ODE shown below

Then a function is created as below for solving the ODE

In the function, following substitutions are made

Function:

function [A] = func_ode(b,g,l,m,t,theta)
theta1 = theta(1);
theta2 = theta(2);
dtheta_dt = theta2;
dtheta1_dt = -(b/m)*theta2-(g/l)*sin(theta1);
A = [dtheta_dt; dtheta1_dt];
end

Coming back to the program, after giving the inputs, angle and time span are entered.

After that, the function is called to solve the ode. The output values are stored in the form of an array ([time,result])

A graph is plotted to visualize the output values.

Further, for loop is created to obtain frames for each output value.

Initial coordinates and final coordinates of the string are defined.

Plotting various plots such as line to create base, string and bob of the pendulum in a single frame. To create a bob, ‘marker’ command is used.

Next step is to create a movie from the frames obtained. For this, ‘movie’ command is used to create a movie.

Code:

clear all
close all
clc

%program to simulate the motion of a pendulum

%Inputs (A per given)
l = 1;  %length in meter
m = 1;  %mass in kg
b = 0.05;    %damping coefficient
g = 9.81;   %Specific gravity in kg/m^3

%initial condition
theta0 = [0;3];

%time points
t = linspace(0,20,500);

%calling function
[time,result] = ode45(@(t,theta) func_ode(b,g,l,m,t,theta),t,theta0);

%plotting graph of the ode
figure(1);
plot(time,result);
xlabel('Time (sec)');
ylabel('Values');

ct = 1;

for i = 1:length(result);
    
    %Initial coordinates
    x0 = 0;
    y0 = 0.5;
    
    %angle of string
    theta_s = result(i);
    
    %End coordinates
    x1 = 1*sin(theta_s);
    y1 = -1*cos(theta_s);
    
    %plotting the values
    figure(2);
    %base on which pendulum is attached
    plot([-0.5 0.5],[0.5 0.5],'linewidth',10,'color','black')
    hold on
    %string of the pendulum
    line([x0 x1],[y0 y1],'linewidth',2,'color','red')
    %bob of the pendulum
    plot(x1,y1,'marker','o','markersize',20,'markerfacecolor','cyan')
    axis([-1.5 1.5 -1.5 1.5]);
    pause(0.03);
    title('Motion of Pendulum');
    hold off;
    M(ct) = getframe(gcf);
    ct = ct+1;
end

%creating movie 
movie(M);
video = VideoWriter('Motion_of_Pendulum','Uncompressed AVI');
open(video);
writeVideo(video,M);
close(video);

%end of program

Outputs:

Figure no. 1

Figure no. 2

Link for the animation : https://youtu.be/KzwozTyb09o

Problems faced:

While typing the program, doubt related to the formula to calculate end coordinates of the string is arrised.

Solved the doubt by searching on the google realted to motion of a pendulum.

No other errors are occured.