Week 3 - Solving second order ODEs

AIM

To write a program in MATLAB to solve the ODE that represents the equation of motion of a simple pendulum with damping.

GOVERNING EQUATIONS

The ordinary differential equation that represents the motion of a simple pendulum is as follows,

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

Where, `theta` is the angular displacement of the pendulum

‘b’ is the damping coefficient

‘m’ is the mass of the pendulum bob

‘g’ is the acceleration due to gravity

‘l’ is the length of the string

Formula to find change in coordinates of the pendulum,

x0 = 0

Y0 = 0 , where x0 and y0 are the fixed point of the pendulum string

`x1 = l*sin(theta)`

`y1 = -l*cos(theta)` , where x1 and y1 are the coordinates that change with respect to change in time.

OBJECTIVE

To create an animation of the pendulum motion.

PROGRAM

All the constants are defined at the beginning of the program , which includes values of b,g,l and m.

Next step is to provide the initial conditions, I.e the initial values of displacement and velocity in an array.

Then we have to specify the time interval.

The syntax for ODE45 for this program,

[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta)

Here ‘t’ is the time , which changes from 0 to 20s.

‘Results’gives us an array with 2 columns containing velocity and displacement values with respect to time.

‘ode45’is the built in MATLAB ODE Solver.

‘ode_func’ is the function program name and its variables given inside the bracket.

‘t_span’specifies the limits for the integration operation. In this case it is 0 to 20.

‘theta’specifies the initial conditions , in this case initial conditions for displacement is 0 and angular velocity is 3 rad/s.

Now we create a for loop to get an animation showing the motion of the pendulum.

Here the animation is created by using the position of the pendulum as a function of time.

After that by using the “movie” command we create an animation .

clear all
close all 
clc

%defining constants values
b=0.05;
g=9.81;
l=1;
m=1;

% specifying initial conditions
theta_0 = [0;3];

% specifying the limits for integration
t_span = linspace(0,20,500);

%ode45 solver 
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

% graph showing velocity and displacement variation w.r.t time
figure(1)
plot(t,results(:,1),'linewidth',3,'color','r')
hold on
plot(t,results(:,2),'linewidth',3,'color','g')
legend('displacement','velocity')
title('graph')
xlabel('time')
ylabel('plot')

%for loop to create an animation of pendulum motion
ct=1;
for i=1:length(t_span)
    theta = results(:,1);
 
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(theta(i));
    y1 = -l*cos(theta(i));
    
    figure(2)
    plot([x0 x1],[y0 y1],'linewidth',3,'marker','o')
    axis([-1 1 -1.2 0])
    title('Motion of the pendulum')
    pause(0.1);
     M(ct)=getframe(gcf);
      ct=ct+1;
    hold on 
    hold off
end

% creating the movie file 
    movie(M)
videofile = VideoWriter('ODE_main_pgm.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

FUNCTION PROGRAM

The function program in this program is named as ‘ode_func’.This contains the equations for representing the simple pendulum motion.

For writing this function program we need to rearrange the ODE mentioned above as follows,

Let `theta = theta_1`

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

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

From this we can find `(d theta_2)/(dt)` as ,

`(d theta_2)/(dt)+(b/m)*theta_2+(g/l)sin theta_1 =0`

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

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

% defining variables
theta1 = theta(1);
theta2 = theta(2);

%equations to be solved
dtheta1_dt = theta2;
dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];

PENDULUM ANIMATION

https://drive.google.com/file/d/1GLBpc3usenLl3cc3y7q1SDBC_UzCdEjB/view?usp=drivesdk

CONCLUSION

The animation showing the motion of the pendulum has been successfully created using the MATLAB program.