Week 3 - Solving second order ODEs

CODE FUNCTION OF ODE 45 FOR THE MAIN PROGRAMME

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

MAIN PROGRAM SOLVING SECOND ORDER ODEs

clear all
close all
clc

% Solving the Second order ODE and show the animation for the second order
% of simple pendulum with damping
% damping coefficient
b = 0.05;
% length of pendulum in m
l=1;
% mass of bob in kg
m = 1; 
%gravity in m/s^2
g=9.81;
% Initial Condition 
theta_0=[0;3];

%Time points

t_span = linspace(0,20,100);

%Solve ODE

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

figure(1)

plot(t,results(:,1)) 
hold on
plot(t,results(:,2))
xlabel('Time')
ylabel('velocity')
legend('Displacement', 'Angular velocity')
grid on 
hold off

displacement=results (:,1); 
velocity=results (:,2);
%fixed coordinates
x0=1;
y0=1;
ct=1;
for i=1:length(displacement)

%coordinates

x1=1-l*sin(displacement(i)); 
y1=1-l*cos(displacement(1));

figure(2)

%plot fixed link
plot([0 2], [1 1], 'linewidth', 10, 'color', 'r')
hold on
plot([1 x1],[1 y1],[x1 x1],[y1 y1],'-O','markersize',15, 'markerfacecolor', 'b', 'linewidth',3)
plot(x0,y0,'-O', 'markersize',3, 'markerfacecolor','b') 
%axis setting
axis([0 2 -0.5 2])
% to change the frame of the movie
M(ct)=getframe(gcf);
ct=ct+l;
hold off
end
%for creating animation video
movie(M)
videofile=VideoWriter('simplependulum.avi','uncompressed AVI');
%open the video file for running the code
open(videofile)
%for writing the vide file
writeVideo(videofile,M)
%for closing the videofile after code
close(videofile)

CODE EXPLANATION

• The given details are assigned as variables in the code
• Initial displacement =0,nitial angular velocity=3
  theta_0=[0;3];
• Then time span of the simple pendulum is set by using linspace command
  t_span = linspace(0,20,500);
• After this we use ode45 function and corresponding 'ode_function' function call in main program
  [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_function(t,theta,b,g,l,m) call ode function
• Then used plot command to plot graph corresponding to time and velocity
• legend command is used for identify each values in the graph.
• For getting the accurate simulation of the simple pendulum corresponding to the second order ODE with damping we use the for loop to create this
• The pendulum has fixed coordinates mentioned as x0,y0 and xl, yl are the other coordinates
• Marker command is used for making bob of the pendulum
• The axis is fixed using axis command
• The counter command is used to make the value changes in each for loop
• Creating an avi file to store the movie and give file name and profile
  videofile=VideoWriter('simplependulum.avi','uncompressed AVI');
• To open the video file for running the code
  open(videofile)
• For writing the vide file
  writeVideo(videofile,M)
• For closing the videofile after code
  close(videofile)

ERRORS