Week 3 - Solving second order ODEs

%second order ode 

clear all
close all
clc

% b = damping coefficient
% l = length of pendulam
% g = gravitational force
% m =- mass

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


% initial condition

%initial displacement = o
%^initial velocity = 5

theta_0 = [0;3] ;


%time points 
%time array for intrgrating ode solution

t_span = linspace(0,20,100) ;

%solve ode

% we have to call a ode function
% ode45 is a inbuilt function
% main inputs are time (t) and angle (theta )
% and so ode is a function of t and theta
% ode_func mean defining a function of input parameters ( t,theta,b,g,l,m)
%calling the function  ode_func( t,theta,b,g,l,m )
% followed by a time array(t_span) for integration  and initial condition (theta_0)



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


% some spelling mistakes are occured in this area
% then solved that errors

plot(t,results(:,1),"Color",'r') 
hold on
% use hold on to show both plots
% otherwise it shows only last plot 
% used different colors for each plot
%labeled that on graph
plot (t,results(:,2),"Color",'b')
xlabel('time')
ylabel ('parameters')
legend('displacement','velocity')

% create a ode_function


  % 2nd order ODE

  % d^2theta/dt^2

  % theta = theta1

  % d(theta)/dt = d(theta1)/dt1 = theta2

  % d^2(theta)/dt^2 = d^2(theta1)/dt^2 = d/dt ( d(theta)/dt )= d(theta2)/dt

  % d^2theta/dt^2 + b/m d(theta)/dt + (g sin(theta)) / l = 0

  % d^2(theta1)/dt^2 + b/m d(theta1)/dt + (g sin(theta12)) / l = 0

  %  d(theta)/dt = (theta2)  ,   d^2(theta1)/dt^2 =  d(theta2)/d

  % d(theta2)/dt + b/m (theta2) + (g sin(theta1)) / l = 0

  % d(theta2)/dt =  - b/m (theta2) - (g sin(theta1))

  % d/dt [ theta1 ; theta2 ] = [ theta2 ; - (b/m )*theta2- (g/l)* sin(theta1)

  % d/dt [ theta1 ; theta2 ] = dtheta_dt




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

used this information to animate the motion of the pendulum

% animation for  the second order of a simple pendulum with damping

% Clear all variables from mengry
% close all charts and figures
% clear command windows

clear all 
close all
clc


%inputs

b = 0.05 ; % b = damping coefficient , 
g = 9.81 ; % g = gravitational force   , m/s^2
l = 1 ; % l = length of pendulam , m
m = 1 ; % m =- mass , kg
 
% initial condition at time=o

theta_0 = [0;3] ; %initial displacement = o , initial angular velocity = 3


%time points 

t_span = linspace(0,20,100) ;  %time array for intrgrating ode solution

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

% ode45(@(t,theta) = find out Angular Displacement with respect to time
% ode_func(t,theta,b,g,l,m) = call ode function
% t_span time array total 100 from 0 to 20


ct = 1;

% initial value for the plot number; set it 1

for i=1:length(results(:,1))

% i is the range of first column of results

x0 = 0; % fixed point of link in x-axis

y0 = -0.05; % Fixed point of link in y-axis

x1 = l*sin(results(i,1)); % angular movement in x-axix


y1 = -l*cos(results(i,1)); % angular movement in y-axix



x2 = (l/2)*sin(results(i,1));


y2= -(l/2)*cos(results(i,1));


% Plot pendulum support line

plot([-0.6 0.6],[0 0], 'linewidth',8, 'color', 'k');

hold on

% hold on use to put next plot in same chart

% Plot pendulum support dot

plot (0,-0.03,'o','MarkerSize', 12, 'MarkerEdgecolor', 'k', 'MarkerFacecolor', 'k');

% Provid Link with length= 1

line([x0 x1],[y0 y1], 'linewidth',1,'color','b')

% Name the link with L=1 by usinf text command
 
t = text((x2+0.04),y2, 'L=1');

s = t.FontSize; % give font size

t.FontSize = 9; % Mention font size of 9

hold on;

% hold on use to put next plot in same chart

% plot pendulum circle

% with markers give size to pendulum

% with markerfacecolor give face color to the pendulum circle 

plot (x1,y1,'-o','markers', 20, 'markerfacecolor', 'r') 

title('Pendulum Damped Animation')

xlabel('x-axis') 

ylabel('y-axis')

grid on

axis ([-1.5 1.5 -1.5 0.5])  % fix axis

pause(0.03) % Pause interval between each frame

hold off % hold off use for turn of the hold on command


% in M(ct), ct represents each and every frame to store in M

% use getframe command to get current frame window in each loop)

M(ct) = getframe(gcf);

% increment of number of frame in each loop

ct=ct+1;

end

% after the end of  for loop , we use a Movie command to geneate a movie realte to
% above program
movie(M);

%create an avi file to store the movie
% we have to give file name and profile 

videofile = VideoWriter('simple pendulam 5.avi' , 'Uncompressed AVI') ;
open(videofile)
writeVideo(videofile,M)
close(videofile)

errors faced in these steps relate to values and spelling. after correcting each error, I got the final output.

youtube link : https://youtu.be/9t_Jryb8LRQ