Week 3 - Solving second order ODEs

Aim:

To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

The ODE function is

Also, to create an animation of the pendulum motion.

Inputs:

b=0.05 (Damping coefficient)

g=9.81 (gravity in m/s^2)

l=1 (Length of pendulum in meter)

m=1 (Mass in Kg)

Initial conditions,

theta_0=[0;5];

Angular Dislplacement is 0 and Angular Velocity is 5

and

t_span=linspace(0,10,500);

Time span is 0 to 10 sec, with 500 intervals.

Solving ODE:

In Matlab ode45 is the inbuilt solver for ODE,

In this, ODE is the function of t and theta.

[t,result]=ode45(@(t,theta) ode_function(b,g,l,m,t,theta),t_span,theta_0)

In ode45 synatx the value of b,g,l,m are taken from the inputs and then we define time span t_span and initial conditions by theta_0

We have called function of name ode_function into it and give a name to ode_function as [dtheta_dt]

First Simplify,

Assume,

theta=theta1 ------>(1)

So, dtheta_dt=dtheta1_dt

dtheta1_dt=theta2 --------->(2)

Then, d^2theta_dt^2=d^2theta1_dt^2

=d_dt(dtheta1_dt)

=d_dt(theta2)

d^2theta_dt^2=dtheta2_dt --------->(3)

So, dtheta2_dt+(b/m)*theta2+(g/l)*sin(theta1)=0 --------->(4)

dtheta1_dt=theta2 --------->(5)

dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1) --------->(6)

Now we have two equations, lets club into matrix form.

d_dt[theta1; theta2]=[theta2; -(b/m)*theta2-(g/L)*sin(theta1)]

and dtheta_dt=[dtheta1_dt; dtheta2_dt]

To find dtheta_dt matrix, we will created the ODE function.

ODE 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

Main code:

clear all
close all
clc

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

%initial condition
theta_0=[0;5];
%time points
t_span=linspace(0,10,500);
%solving ODE
[t,result]=ode45(@(t,theta) ode_func(b,g,l,m,t,theta),t_span,theta_0);

plot(t,result(:,1))
hold on
plot(t,result(:,2))
ylabel('displacement and velocity')
xlabel('time')
grid on
legend('displacement','velocity')

%to create animation
ct=1;
for i=1:length(t)
    x0=0;
    y0=0;
    
    x1= l*sin(result(i,1));
    y1=-l*cos(result(i,1));
    plot([-2 2],[0 0],'linewidth',3,'color','r');

    line([x0 x1],[y0 y1],'linewidth',4);
    
     axis([-3 3 -3 3]);
     
     hold on
     plot(x1,y1,'o','markers',20,'markerfacecolor','y');
     hold off
     
     grid on
     pause(0.03);
     
     
     M(ct)= getframe(gcf);
     
     ct=ct+1;
end

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

M is used to capture each frame when the values of the plottings are changing since for is used.

This value of M is further used in creating the video/movie.

A video file is created and then these M values are loaded in the video file through the command 'writeVideo'.

Output:

The below plot shows the change in the velocity and the displacement of the pendulum.

The video showing the output of the second-order ODE function is shown in the below youtube video link:

Link:- https://youtu.be/PwDot2NX2Cs