Solving Second Order ODE's

AIM- To write a Matlab program to solve the following ODE and simulate pendulum motion.

`(d^2theta )/dt^2+ b/m*(d theta)/dt+g/(l*sin theta)=0`

fig 1 schematic diagram of a simple pendulum

Theory-

A simple pendulum is the one which works on newtons second law of motion , which consists of a small ball hanging to a mass less string of length 'L' fixed at a point P. When displaced to an initial angle and released the pendulum will swing back and forth with periodic motion.

The below equation represents the ODE of a pendulum.

`(d^2theta )/dt^2+ b/m*(d theta)/dt+g/(l*sin theta)=0`

`b=` damping co-efficient

`m=` mass of ball in kg

`g=` gravity in m/s^2

`l=` length of pendulum in m

Since the above equation is of second order differential equation we need to break into two first order differential equations and solve as follows-

`(d^2theta )/dt^2+ b/m*(d theta)/dt+g/(l*sin theta)=0 (a)`

let `theta=theta_1`

`(d theta)/dt=(d theta_1)/(dt)=theta_2`

Now diff. w.r.t t

`(d^2 theta)/dt^2=(d^2 theta_1)/(dt^2)=(d theta_1)/dt` sub in equation 'a'

Now we have two first order differential equations

`(d theta_1)/(dt)=theta_2` -eq '1'

by re arranging eq 'a' we get

`(d^2theta )/dt^2=- b/m*(d theta)/dt-g/(l*sin theta)` - eq '2'

these equations can also be represented as

`(d)/(dt)((theta_1),(theta_2)) = ((theta_2),(- b/m*(d theta)/dt-g/(l*sin theta)))`

Software used- MatlabR2018a

INPUT MAIN PROGRAM-

clear all
close all 
clc

b=0.05;
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_function(t,theta,b,g,l,m),t_span,theta_0);
plot(t,results(:,1))
hold on
figure (1)
plot(t,results(:,2))
grid on
xlabel('time')
ylabel('plot')

% animation
ct=1;
for i=1: length(results)
    theta=results(i);
    
    x0 = 0;
    y0 = 0;
    x1 = l*sin(theta);
    y1 = -l*cos(theta);
    figure(2)
    plot([-1 1],[0 0], 'linewidth',4,'color','r');
    axis([-2 2 -2 2]);
    line([x0 x1],[y0 y1],'linewidth',2,'color','g');
    hold on;
    plot(x1,y1,'o','markersize',20,'markerfacecolor','b','markeredgecolor','b');
    grid on
    hold off;
    M(ct) = getframe(gcf);
    ct=ct+1;
    

end

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

Matlab code explaination-

//

clear all
close all 
clc

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

Initially the values of , damping co-efficient (b)=0.05,gravity(g)=9.8 m/s^2 , length of pendilum(l)=1 meter, mass of bob(m)=1kg are given .

//

% initial conditions
theta_0=[0,5];
% time points 
t_span=linspace(0,10,500);
% solve ode
[t,results]=ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,theta_0);

Now the initial conditions are determined and stored in an array theta_0 , initial displacement is zero and initial velocity is 5m/s^2.

Next t_span a time array is created and given time from 0 to 10 seconds , of which 500 values between them are generated using linspace command.

Further an inbuilt solver ode45 which solves 1st order diffrential equations is assigned and solved w.r.t t and theta as an anonymous function using @(t,theta),the nature of ode45 solver is that it accepts only an indipendent and dependent variable ,but we used @(t,theta) to act as an anonymous function so that it enters into the function ptogram ode_function which has different variables and solves it.

//

Now ode45 solver enters the function program

INPUT FUNCTION CODE-

function [dtheta_dt]= ode_function(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

If theta is differentiated once we get displacement and it is diffrentiated twice we get velocity.

Further theta1 is assigned to theta(1) and theta2 is assigned to theta(2).

Now to solve displacement i.e 1st ode, dtheta1_dt=theta2 equation is evaluated and to solve velocity dtheta_dt =-(b/m)*theta2-(g/l)*sin(theta1) is evaluated.

Then evaluated [dtheta1_dt;dtheta2_dt] is stored in an array dtheta_dt seperately in two coloumns.

Now the function code get ended, and goes to main program and stores the values in [ t, results ].

//

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

Now a plot is generated displacements and velocities values using plot command under the name figure(1).

Further x axis and y axis labels are given to the plot usng xlabel and yabel commands.

//

Consider the initial co-ordinates rigid support of the pendilum as(0,0) and if we resolve we get x=l*sin(theta) as x co-ordinate of end poin and y co-ordinate as y=-l*cos(theta) since it is in -y axis.

Now a for loop is generated for generating animation.

ct=1;
for i=1: length(results)
    theta=results(i);
    
    x0 = 0;
    y0 = 0;
    x1 = l*sin(theta);
    y1 = -l*cos(theta);
    figure(2)
    plot([-1 1],[0 0], 'linewidth',4,'color','r');
    axis([-2 2 -2 2]);
    line([x0 x1],[y0 y1],'linewidth',2,'color','g');
    hold on;
    plot(x1,y1,'o','markersize',20,'markerfacecolor','b','markeredgecolor','b');
    grid on
    hold off;
    M(ct) = getframe(gcf);
    ct=ct+1;
    

end

for loop is exicuted for i=1 to obtained theta values i.e displacement values.

since the one sine of the pendilum has fixed end it is represented as x0 = 0;y0 = 0; then the oscillating end is represented as x1=l*sin(theta),y1=-l*cos(theta), for every value of change in angle i.e theta we get a new position, to store these positions in frames created an array M(ct) ,getframe(gcf) command stores all the frames in the array M.

A variable ct is used to store the positions of each frame , ct value iterates untill the final frame is stored, for which we compile all the frames and create an video file and run it.

To represent the fixed end a red coloured line with linewidth 4 is obtained using plot command.

Also a line coloured in green is joined from fined point(0,0) to the oscillating end using the line command.

To create a bob at the end we use plot command and marker size as 20 and color as blue.

axis is set from -2 to +2 in both x and y axis.

//

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

using the command VideoWriter all the frame gets compiles into a video and stores as pendulam_ode.avi in an array videofile.

Now we open the videofile using open(videofile);

Then play the video using writeVideo(videofile,M) command.

Then close it using the command close(videofile).

//

OUTPUT-

The displacement and velocity of the pendulum are obtained at different time intervels and the animation is generated to examinethe motion sucessfully.

fig 2 represents displacement and velocity at diffrent time interevels

fig 3 represents a frame of pendulum simulation

Animation-

The animation can be viewed using the below mentioned link.

https://youtu.be/1wbvm_jmJLs

Conclusion-

Hence, the second order ODE is solved and simulation of simple pendulum was created sucessfully.