Solving second order ODE using MATLAB

Aim: To generate an animation of a simple pendulum in motion using MATLAB.

Introduction: Let us assume a spherical bob is suspended with a string and the dimensions are as follows-

Length of pendulum string from support to bob = L = 1m.

weight of bob = m = 1kg.

damping cooefficient = b = 0.05

gravity = 9.81`m/(sec)^2`

The differential equation for simple pendulum is given as follows:

`(d^2theta)/dt + b/m* [d theta]/dt`` + g/L*sintheta =` 0

Converting second order into 2 first order ODE:

Suppose:

`theta = theta1` so, `[d theta_1]/dt = theta_2` ===> consider this as #1`

`therefore` `(d^2theta)/dt = (d^2theta_1)/dt^2 = d/dt([d theta_1]/dt) = [d theta_2]/dt` ===> consider this as #2`

`(d^2theta)/dt + b/m* [d theta]/dt + g/L*sintheta = 0`

`therefore [d^2theta_1]/dt + b/m*(d theta_1)/dt + g/L*sintheta_1 = 0`

`therefore (d theta_2)/dt + b/m * theta_2 + g/L*sintheta_1 = 0`

`therefore (d theta_1)/dt = theta_2 and (d theta_2)/dt = -b/mtheta_2 - g/Lsintheta_1`

Now representing this in matrix format, we get

`d/dt[[theta_1],[theta_2]] = [[theta_2],[-b/mtheta_2-g/Lsintheta_1]]`

`therefore Theta =``[[theta_2],[-b/mtheta_2-g/Lsintheta_1]]`

Coding in MATLAB:

Now the value `Theta` is saved in different file as "ode_func" function and we will call this function when we want while coding for the project.

In the code, then we specify the dimensions of the bob,string and values of associated parameters

then we specify the initial condition, time span by using "linspace" command

"ode45" solver is used to solve the higher order odes.

ODE45 EXPLAINED

ode45 is an in-built solver in MATLAB which takes the input arguement, function definition,time span values,and initial condition and gives us in return time, angular displacement and angular velocity in matrix form.

@(t,theta) = t and theta are input arguements,

ode_func = is the name of the function we saved separately and it defines the function by possessing the values of parameters.

t_span = time span,

theta_0 = initial condition.

PLOTTING:

After we solve the ODE ,we move forward to generate the plots.

"subplot" command is used to fit multiple plots in one frame by calling out their positions.

Plots representing angular displacement and angular velocity are plotted.

And a GIF for representing the motion of simple pendulum is generated. At first "for loop" and "getframe" command is used to produce plots at different values and then they are stiched together.

"Movie" command is used to generate a movie by calling the coded .m file.

Coding file is as follows:

clc
close all
clear all

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

% 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);

%plot(t,results(:,1))
%hold on
%plot(t,results(:,2))
%ylabel('plot')
%xlabel('time')


%% code for animation

subplot(5,4,[13,14,17,18]);
plot(t,results(:,1),'color','b');
hold on;
plot(t,results(:,2),'color','r');
hold off;
subplot(5,4,[15,16,19,20]);
plot(results(:,1),results(:,2));
xlabel('position-radians');
ylabel('angular velocity');


ct=1;
for i = 1:length(results(:,1));
    x0=0;
    y0=0;
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);
    plot([-1 1],[0 0],'linewidth',4,'color','k');
    hold on
    line([x0 x1],[y0 y1],'linewidth',2,'color','b');
    hold on;
    
    
    plot(x1,y1,'-o','markers',20,'markerfacecolor','r');
    
    grid on
    axis([-1.5 1.5 -1.5 0.5]);
    pause(0.0003);
    hold off
    M(ct)=getframe(gcf);
    %Text=sprintf("pendulum%d.png",ct);  to give name in octave
    %saveas(1,Text);   saveas is command (octave)
    ct=ct+1;
    
end

movie(M)
videofile = VideoWriter('ode_main.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

the ode function code is as follows:

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

workspace for the code is as follows:

in the results file we get values for angular displacements and angular velocity.

Animation file for the simple pendulum is follows:

In order to dampen the oscillation, the damping coefficient can be increased and the animation video then is as shown bellow.

damping coefficient = b = 0.1

ERRORS while coding:

the suspended pendulum is supposed to be hanging down. so the plots has to be in third and fourth quadrant.

Solution:

so in order to change the quadrant "minus" sign had to be added for the y-component of the string line.

Hence, the Simple Pendulum is represented using ode45 solver and Movie commands in MATLAB.