Week 3 - Solving second order ODEs

Aim :-

To solve a transient behavior of a system. (Simple Pendulum).

Condition to consider :-

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Output:-

Output Video:-

Code:-

% Solving differential equations

clear all
close all
clc

% Inputs

b = 0.05; % Damping coefficient
g = 9.81; % gravity
l = 1; % length
m = 1; % Mass

% Initial condition
theta_0 = [0;3];
    
% Time points
t_span = linspace(0,20,500);

% Solve ODE
[t, results] = ode45(@(t,theta) ode_func(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('Angular velocity & Angular Displacement')
grid on
legend('Angular Displacement','Angular Velocity')
title('Variation of Angular displacement over Angular velocity with respect to time')

% Simple pendulum animation
theta2 = results(:,1);

ct = 1; % Counter for looping

for i = 1:length(theta2);
    THETA2=theta2(i);
    
    x0 = 0; % Initial Position along X axis
    y0 = 0; % Initial Position along Y axis
    
    x1 = l*sind(THETA2); % Position w.r.t to theta
    y1 = -l*cosd(THETA2); % Position w.r.t to theta
    
    % Plotting
    figure(2)
    plot([-0.05, 0.05], [0,0],'color','b','linewidth',6)
    hold on
    plot([x0, x1], [y0, y1],'color','r','linewidth',3)
    hold on
    plot(x1,y1,'o','markers',25,'markerfacecolor','y')
    hold off
    grid on
    title('Simple Pendulum Animation')
    M(ct) = getframe(gcf);
    ct = ct+1;
end

% Movie making 
movie(M)
videofile = VideoWriter('simple_pendulum','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function code to solve ODE

% Function ODE

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