Week 3 - Solving second order ODEs

function [dtheta_dt] =ode_fun(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

clear all
close all
clc

b = 0.05; 
g = 9.81;
l = 1;
m = 1;
%initial condition
theta_0 = [0;3];
%time point
t_span = linspace(0,20,500);
%solve ODE
[t,results] = ode45 ( @(t,theta) ode_fun(t , theta , b , g , l , m),t_span,theta_0);
ct =1;
for i= 1:length(t)
         x0 = 0;
         y0 = 0;

         x1 = l*sin(results(i,1));
         y1 = -l*cos(results(i,1));   
         figure(1)
         plot([-1 1],[0 0],'LineWidth',4,'color','y');
         axis([-2 2 -2 2]);
         line([x0 x1],[y0 y1],'linewidth',1,'color','b');
         hold on,   
         plot(x1,y1,'O','markersize',20,'MarkerFaceColor','r');
         grid on;
         hold off,
  
         M(ct) = getframe(gcf);
         ct=ct+1;

end
movie(M)
videofile = VideoWriter('Simple_pendulam','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

clear all, close all, clc command is used for clearing the workspace, command window, and for closing all graphical tabs.

STEP1: Assigning the input values and ie. damping coefficient, gravity, length of the pendulum, and mass of the pendulum (b = 0.05,g = 9.81,l = 1,m = 1) as given in the problem.

STEP2: Assigning the initial conditions for the ODE (displacement and velocity) into a variable 'theta_0' ie theta_0 = [0;3] as provided in the problem. For integrating the function time array is required which can be assigned to the variable 't_span', As we need multiple values for the array linspace command is used, The time range 0-20 is given in the problem and we assign 500 values (t_span = linspace(0,20,500)).

STEP3: To solve the ODE we will be needing two functions , the solver and the other to convert the 2nd order ODE into a set of first-order ODEs. For the first function, we have an inbuild Matlab function ie 'ode45' and for the second we will assign the theta values and equation in variables and compile them into a single array ie. 'dtheta_dt = [dtheta1_dt; dtheta2_dt] where dtheta1_dt = theta2,dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1) and theta1 = theta(1) , theta2 = theta(2)'.The function name and inputs assigned will be 'ode_fun(t , theta , b , g , l , m)'.

STEP4: solving the ODE,assigning the dependend and independend variable of the ODE ie. t and theta into the solver ode45 then calling the function 'ode_fun(t , theta , b , g , l , m)' assigned with input parameters after that we provide the integration range (t_span = linspace(0,20,500))array followed with initial condition 'theta_0'.Assigning all this to a variable array 't'and 'results' ie [t,results] = ode45 ( @(t,theta) ode_fun(t , theta , b , g , l , m),t_span,theta_0)

STEP5:As we have solved the ODE, now we have to present it in a graphical form with the help of pendulum motion, and to do so we have to take the frame by frame pendulum position by assigning a counter (ct=1) and loop.' for' command is used for the looping process where 'i' is assigned with the number of integration 't' ie 1:length(t) then assigning the initial x0 = 0 and y0 = 0 points after that assigning the final point x1 and y1 which can be done by multiplying the length of the pendulum with angular results displacement (theta values) ie x1 = l*sin(results(i,1)) similarly y1 = -l*cos(results(i,1)) but as we need hanging oscillation so we use negative length. The rest is to represent the supporting rod and length of the pendulum with the help of 'plot' and 'line' command ie. plot([-1 1],[0 0],'LineWidth',4,'color','y'),line([x0 x1],[y0 y1],'linewidth',1,'color','b') we use 'hold on' command to hold the pendulum length frame so that bod can be assigned later . then we represent the bob of the pendulum with the help of markersize and marker facecolor command in plot command which creates the sphere ie. plot(x1,y1,'O','markersize',20,'MarkerFaceColor','r') then we use the hold off command so that we can have a smooth transition without glitch effect of the pendulam length. then with the help of M(ct) = getframe(gcf) we can store the frame after that increment the counter(ct=ct+1) so that we can get frames at the end of every value of displacement,close the loop with the help of end command.

STEP6:creating the animation, With the help of movie command we assign the frames stored in 'M' ie 'movie(M)'.then using videowriter naming and file formatting by assigning it to a variable videofile (videofile = VideoWriter ('Simple_pendulam' ,'Uncompressed AVI')) next process is to open the file and with the help of 'writevideo' command assign the frames into it then close the file (open(videofile),writeVideo(videofile,M),close(videofile).

OUTPUT

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

the animation show the pendulum oscillation

PROBLEM FACED.

I had assigned the displacement values in sind() command and I was not able to get the proper oscillation, later while reassuring the values I found out that all values are in radian then I used sin() command.

I didn't assign the negative length for y1 and I was facing the upside-down pendulum oscillation I tried adding negative inside the results but it was giving an error that cant assign a negative value. after that I tried adding it to the length and it worked.