Week 3 - Solving second order ODEs

Main code:

clear all
close all
clc


% Inputs 
b= 0.05; % damping coefficient
g= 9.81; % acceleration due to gravity (in m/s^2)
l= 1; % length of the spring (in m)
m= 1; % mass of the bob (in kg)

% initial conditions
% angular displacement=0, angular velocity= 3 rad/s
theta_0= [0;3];


% timepoints
% Oscillation time between 0-20 secs
t_span= linspace (0,20,300);

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

%plotting
figure(1)
hold on
plot (t, results (:,1))
plot (t, results (:,2))
title('Damped motion of a pendulum');
legend('Angular displacement', 'Angular velocity');
ylabel ('plot')
xlabel('time')
 

% for loop for animation 
ct=1;

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

% fixed end
x1= 0;
y1= 0;

% position of bob
x2=  l*sin(results(i,1));
y2= -l*cos(results(i,1));
figure(2)
plot([-1 1],[0 0],'color','r','linewidth',3.5);
hold on
plot([x1 x2],[y1 y2],'linewidth',2.5);
hold on
plot(x2,y2,'o',Markersize=30,Markerfacecolor='black',Markeredgecolor='black');
axis([-1.5 1.5 -1.5 1]);
title('Simple pendulum motion');
pause(0.03);
hold off

M(ct)=getframe(gcf);
ct=ct+1;

end

% Creating animation
movie(M)
videofile=VideoWriter('Pendulum animation.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function code:

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

Explanation:

First step is to define the input datas; damping coeff., acc. due to gravity, length of the string and mass of the bob. Secondly, is to define initial conditions, angular velocity, angular displacement along with the timespan. Here, we are assuming angular displacement= 0 and angular velocity= 3. The function code is used to perform a speciifc task. Here it is used to solve ODE. It is in the format of .m. The function definition is declared in the main code itself and its value is taken from the function code. To solve the ODE, ODE 45 is being used. It is a function that solves differential equations. Here it is calling 3 variables, t_span, [t,t_results] these will be integrated and calculated w.r.to time and plotted in the form of a graph. The graph will be angular displacement vs time and angular velocity vs time. The curves indicating the respective parameters are indicated with different colors as in the graph.

Graph:

The second part of the program talks about the animation of the bob. Here for loop is used. Its purpose is to execute a particular intended work number of times. Here timespan values are assigned to variable 'i'. Here initial conditions of the bob are assumed as x1=0, y1=0 & x2, y2 are calculated using the trignomteric functions. The purpose of marker command is used to create the bob shape of the pendulum. Pause command is used to pause the command for certain time i.e 0.03s. M command is used to capture the frame every time when for loop runs. Here ct variable is used, each time the for loop runs it is incremented by 1. Finally, it is ended using End command. Finally, animation is created. VideoWriter command is used to assign the video format. Here .avi is the format used. Open command is used to open the video. writeVideo command is used to write the video. Finally animation is closed using close command.

Outputs:

Animation link:

https://youtu.be/CEek9zYLaBI

Error while working on the program:

1. As matlab is case sensitive, VideoWriter command has 'V' , 'W' as upper cases

Remedy:

After giving the command in matlab, try to check for syntax. In this way we can avoid this error