Week 3 - Solving second order ODEs

Solving second order ODE using MATLAB

Introduction: A second-order derivative is a derivative of the derivative of a function. It is drawn from the first-order derivative. So we first find the derivative of a function and then draw out the derivative of the first derivative.

Aim: To write a program that can solve the given ODE that represents the motion of the simple pendulum with damping and create the animation of the motion

Theory: A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging back and forth. The time for one complete cycle, a left swing and a right swing, is called the period. The period depends on the length of the pendulum and also to a slight degree on the amplitude, the width of the pendulum's swing.

Code for Solving Second order ODE

clear all
close all
clc
 
b=0.2;
g=9.81;
l=1;
m=0.1;
 
%Initial conditions
theta_0=[0;5];
 
%Time points
t_span=linspace(0,10,500);
 
%solve ODE
[t, results]= ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
 
%Plotting
plot(t,results(:,1));
hold on
plot(t,results(:,2));
ylabel('Plot')
xlabel('Time')

ct = 1;
 for i = 1:length(results(:,1))
     x0 = 0;
     y0 = 0;
     x1 = l*sin(results(i,1)+pi);
     y1 = l*cos(results(i,1)+pi);
     
     plot([x0 x1],[y0 y1],'-o','MarkerSize',5,'MarkerFaceColor','g','MarkerEdgeColor','g','linewidth',3)
     axis([-2 2 -2 2])
     hold on
     plot([-1 1],[0 0],'linewidth',3,'Color','[0 0.2 0]')
     M(ct) = getframe(gcf)
     hold off
     ct = ct+1;
 end
 
 %Movie
 movie(M)
 v_f = VideoWriter('Simple_Pendulum.avi','Uncompressed AVI');
 open(v_f)
 writeVideo(v_f,M)
 close(v_f)

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 :

• Length, Mass, Damping coefficient and Acceleration due to gravity are given as fixed parameters.
• Initial time and angle of the bob is given.
• Using ode45 command ODE function is introduced.
• Velocity and Displacement is Plotted.
• Using for loop the mass (bob) at different points and time are marked and plotted.

Plot:

Youtube link for the pendulum animation

https://youtu.be/vhG2hBKbNXI

Drive Link for Report

https://drive.google.com/file/d/1ULW1LaR1L20rI8_8fS9BLoF1cYI55rH9/view?usp=sharing