Week 4.1 - Solving second order ODEs

Solving second order ODEs

Aim:

To solves the ODE which represents the equation of motion of a simple pendulum with damping.

Objective:

Write a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping and create an animated video of the output obtained by solving this ODE.

Theory:

The motion of a simple pendulum is a basic classical example of simple harmonic motion, consisting of a small bob and a massless string. Movement of the pendulum depends on Newton's second law. When this law is written down, a second-order Ordinary Differential Equation is formed that describes the position of the end of the pendulum w.r.t time.

Equation of motion is:

In the above equation,

• g = gravity in `m/s^2`.
• L = length of the pendulum in the meter.
• m = mass of the ball in kg.
• b = damping coefficient.

Since the given function is in the second-order differential equation we will convert it into first order to solve the ODE.

Given Data:

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

MAIN CODE:

%second order ODEs

clear all
close all
clc

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

% initial condition
theta_0 =[0;3];

% time 
t_span = linspace(0,20,500)

% slove simple pendulum using ode equation
[t x] = ode45(@(t,theta) pendulum_2nd_order(t,theta,b,l,g,m),t_span,theta_0);

% ploting
subplot(2,1,1)
plot(t,x(:,1),'LineWidth',3);
xlabel('time t in sec');
ylabel('angular displacement');
grid on
hold on
subplot(2,1,2)
plot(t,x(:,2),'LineWidth',3);
xlabel('time t in sec');
ylabel('angular velocity');
grid on

% animation

X0 = 0;
Y0 = 0;

ct = 1
for i =1:length(t)
    X1 =l*sin(x(i,1));
    Y1 =-l*cos(x(i,1));

    figure(2)
    plot ([-1 1],[0 0],'linewidth',6,'color','g');
    axis([-2 2 -2 2]);
    line([X0 X1],[Y0 Y1],'linewidth',2,'color','b');
    hold on;
    plot(X0,Y0,'O','markersize',3,'MarkerFaceColor','y');
    plot(X1,Y1,'O','markersize',30,'MarkerFaceColor','r');
    title('damping of simple pendulum');
    legend('support rod','length of pendulum','fixed point','bod of pendulum');
    grid on; 
    hold off;
    M(ct) =getframe(gcf);
    ct =ct+1;
end

movie(M);
videofile=VideoWriter('simple_pendulum','Uncompressed AVI');
open(videofile);
writeVideo(videofile,M); 
close(videofile);

FUNCTION CODE:

function [dtheta_dt] = pendulum_2nd_order(t,theta,b,l,g,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

• In the first step Input and initial condition of a simple pendulum are fed.
• Then the program is written to solve the ODE in which MATLAB inbuilt function (ode45) is used with constant theta, t, and another function is used, which is defined separately in function program.
• In the next step result of solved ODE is plotted w.r.t time using plot command.
• Now iteration of co-ordinate is fed using for loop to create animation.
• Then plot command is used to plot pendulum string and bob.
• In the next step, the move command is used to create an animated video of the oscillating pendulum.

OUTPUT:

Plot output:

The output obtained is graphical representation velocity curve and displacement curve w.r.t Time.

Angular Displacement and Angular Velocity plots were generated for a time range of 20 sec.

ANIMATION:

Conclusion:

• Using the MATLAB function ordinary differential equations have been created.
• Hence the simple pendulum animation is created.