Challenge

Week 3 - Solving second order ODEs

HHarsha Korlapati14 Jan 2022

AIM - To solve second order ODE to simulate the transient motion of a simple pendullum and create animation file.

MATLAB CODE-

clear all
close all
clc

b = 0.05;  % damping coefficient
g = 9.81;  % gravity in m/s^2
l = 1;     % length of the pendullum in 'm'
m = 1;     % mass of the ball in kg's

%inputs
theta_a = [0;5];

%% TIME SPAN
t_span = linspace(0,20,500);

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

%plot
figure(1)
plot(t,results(:,1))
hold on ;
plot(t,results(:,2))
legend('time vs Displacement','time vs velocity')
xlabel('Plot')
ylabel('Time')


%% CODE FOR ANIMATION AND VIDEO FILE

ct=1;
for i = 1: length(t)
x0 = 0;
y0 = 0;
THETA = results(i,1);
x1 = l*sind(THETA);
y1 = -(l*cosd(THETA));
figure(2)

plot([x0 x1],[y0 y1]);
hold on 
plot(x1,y1,'marker','.','markersize',40)
line([-0.1 0.1],[0,0],'linewidth',2,'color','k')
hold off
axis([-0.2 0.2 -1.5 0.2 ])
 pause(0.03)
 M(ct)= getframe(gcf);
  ct = ct+1;

end 

movie(M)
videofile = VideoWriter('simple_pendullum')
open(videofile)
writeVideo(videofile,M)
close(videofile)

STEP BY STEP PROCEDURE-

To solve the second order ode we have to break it to first order ode.
The first order ODE's are dtheta1_dt and dtheta2_dt
Assign them to the solution vector.
Here we are using ODE45 solver to solve the ode function.
[t,results] = ode45(@(t,theta)ode_fun(t,theta,b,g,l,m),t_span, theta_a)

here t_span = [t0 tf] integrates the system of differential equations y' = f(t,y) from t0 to tf

ode_func is a fuction created damping coefficient , length of pendullum, gravity and mass of the ball and also breaking second order function to first order code also carried out in this function.
Then plot the results theta1 and theta2 against t the we will get the plot of displacement vs time and velocity vs time.
fix x0 and y0 as (0,0).
From basic trignometric formulae we get x1 and y1 values.
Vary x1 and y1 values with respect to theta
Use marker command to create the ball.
finally create the movie accordingly.

FUNCTION-

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

PLOT -

Plot for time vs displacement and time vs velocity