SECOND ORDER ODE Using Matlab

ODE is used to describe the transient behavior of the system. the way pendulm moves depends on the newtons second law

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

To write a program in Matlab\\Octave that will simulate the pendulum motion, just like the one shown in the start of this challenge.

use,

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.The ode45 command is used to integrate the values using the function odefunc by specifing the time interval and initial conditions for the function to start off.

MAIN CODE

clear all
close all
clc

b=0.05;
g=9.81;
l=1;
m=0.1;

%initial conditiong
theta_0= [0;3];

%timepoints
t_span = linspace(0,10,500);

%solve ODE
[t, results]= ode45(@(t,theta) ode_func(b,g,l,m,t,theta),t_span,theta_0);
figure(1)
plot(t,results(:,1),\'color\',\'r\')
grid on
hold on
plot(t,results(:,2),\'color\',\'b\')
xlabel(\'time\')
ylabel(\'angular velocity & angular displacement\')
legend(\'angular Displacement\',\'angular velocity\')
title(\'2nd order ODE\')


%creating movie

ct=1
for i= 1:length(t)
    x1=0;
    y1=0;
    x2= 1*sin(results(i,1));
    y2= -1*cos(results(i,1));
    x3=-1;
    y3=0;
    x4=1;
    y4=0;
    figure(2)
    plot([x3 x4],[y3 y4],\'linewidth\',2,\'color\',\'r\')
    hold on
    plot([x1 x2],[y1 y2],\'linewidth\',2,\'color\',\'r\')
    hold on
    plot(x2,y2,\'o\',\'markersize\',15,\'markerfacecolor\',\'y\')
    axis([-2 2 -2 0])
    hold off
    pause(0.3);
  
    M(ct)=getframe(gcf);
    ct=ct+1;
end

movie(M);
videofile= VideoWriter(\'secondorder_ODE.avi\',\'uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);

The ode45 command is used to integrate the values using the function odefunc by specifing the time interval and initial conditions for the function to start off.

FUNCTION CODE

function [dtheta_dt] = ode_func(b,g,l,m,t,theta)

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

ERROR FACED DURING PROGRAM

AXIS Error

Solved by changing axis to (-2 2 -2 0)

PLOT Angular Velocity Vs Displacement

Youtube link of Animation of Pendulum

https://www.youtube.com/watch?v=P964JWc5Mf0&feature=youtu.be