Week 3 - Solving second order ODEs

Objective : To write a program to execute the behaviour of a simple pendulum & create an animation file.

A pendulum is a weight suspended from a pivot so that it can swing freely. When a pendulum is displaced sideways from its resting equilibrium position, it is subjected to a restoring force due to gravity that will accelerate it back towards the equilibrium position

Consider a pendulum with mass 'm' , length 'l' , making an angle `theta`,when resolved into forces as per mentioned in the equation above we get mgsintheta & mgcostheta`

This ODE below represents the equation of motion of a simple pendulum with damping :

where ,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2

An ordinary differential equation is a differential equation containing one or more functions of one independent variable and the derivatives of those functions.

The below code is the function code , that basically breaks down the second order differential equation to 2 simple first order differential equation's.

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

The next step is to give the inputs of b,g,l,m, theta .

• We give the initial condition of theta_0 , we use 0 as the initial displacement & velocity as 5 rad/sec.
• We use linspace command to get 500 values between 0 & 10 seconds.
• Solve ode is a function, ode45 is a standard function in MATLAB which contains dependent & independent variable's along the function.
• ode_func defines the function & its representation.
• @ Means call the function for the given second order differential equation w.r.t time & initial conditions
• We plot the displacement & velocity graph.
• To start with the animation we use the for loop whose length ranges fom 1 to x1 , y1 position of the pendulum.
• Pendulum is a string plotted using array brackets & ball of pendulum is ploted using marker o.
• The axis is defined & pause command maintains a time interval of 0.01 second between each value of i , Hold off command is used to remove the plot to have only one plot on the screen & hold on is used to draw all elements of the plot.
• Video(M) command is used to make the avi file & write the video after simulation.

close all
clear all
clc
b = 0.05;
g = 9.81;
l = 1;
m = 0.1;

%initialcondition
theta_0 = [0;5];

%times point
t_span = linspace(0,10,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));
xlabel('Position')
ylabel('Angular velocity')
grid on;

videofile = VideoWriter('pendulum_animation.avi','uncompressed AVI');
open(videofile)

%Simulation
ct = 1;
for i = 1 : length(results(:,1))
    x0 = 0;
    y0 = 0;
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    
    plot([-0.5,0.5],[0,0], 'linewidth',3,'color','r');
    line([x0 x1], [y0 y1], 'linewidth',2,'color','g');
    hold on ;
    
    plot(x1,y1, 'o', 'markers',20,'markerfacecolor','b');
    grid on
    axis([-1 1 -1 1])
    pause(0.01);
    hold off
    M = getframe(gcf);
    writeVideo(videofile,M)
    ct = ct + 1;
end

%createing Movie
%movie(M)


close(videofile)

https://youtu.be/pjWCo3H_Z6U

Conclusion : The program & animation file has been executed as per the question.

We did encounter error's but were only typo errors spelling errors to be precise & were solved in weekly doubt sessions in MATLAB meetings