Week 3 - Solving second order ODEs

AIM : To write a Matlab program to solve the ODE which represent the simple pendulum motion and the create animation file of the simple pendulum.

( Using Octave 4.0.0 GUI )

OBJECTIVE OF THE PROJECT:

Our main objective is to solve ODE of the simple pendulum with a damper and to stimulate the motion between 0 to 20 seconds for angular velocity 3 rad/sec, angular displacement 0 at time zero sec and to get plot of time vs (angular displacement and angular velocity), animation of a simple pendulum.

ABOUT SIMPLE PENDULUM :

A pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging it back and forth.

Equation of motion of a simple pendulum with a damper is given by :

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

Where, b = Damping coefficient in N-s/m.

m = Mass of the ball in kg.

l = Length of the pendulum in meter.

g = Acceleration due to gravity m/s².

Let,

`theta=theta1`

differentiating it with respect to time we get

`(d theta)/(dt)=(d theta1) /(dt)=theta2` - (eqn 1)

Again differentiating it with respect to time we get

`(d^2 theta)/(dt^2)= (d^2 theta1)/(d t^2)=(d)/(dt){(d theta1)/(dt)} =(d)/(dt)(theta2)`

We know, the eqaution of motion of simple pendulum with a damper is:

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

This can be written as

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

Again we can write the previous equation as

`(d theta2)/dt+(b)/m theta2+(g)/l sin (theta1) =0` - (eqn 2)

From eqn 1 and eqn 2 we get two ODE :

`(d theta1)/dt=theta2`

and

`(d theta2)/dt=-{(b)/m}theta2-sin(theta1)`

Now writing these ODE in matrix form

```(d)/dt[[theta1],[theta2]]=``[[theta2],[-(b)/mtheta2-(g)/lsintheta1]]`

By using these matrix we can solve the eqn of a simple pendulum

Code for calling the function:

function [dtheta_dt] = damping_ode(theta,t)
   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

Code to solve ODE :

% Solving 2nd order ODE
clear all
close all
clc

%inputs

% Damping coefficient in Ns/m
b= 0.05;
% Acceleration due to gravity in m/s^2
g= 9.81;
% Length of  the string in meter
l= 1;
% Mass of the ball in kg
m= 1; 

%initials conditions 
theta_0 = [0;3];

%time points
t_span = linspace(0,20,350);

%solving ODE
[results] = lsode('damping_func',theta_0,t_span);

%plotting
figure(1)
plot(t_span,results(:,1))
hold on
plot(t_span,results(:,2),'color','r')
xlabel('Time')
ylabel('Angular velocity and Angular displacement')
legend('Angular displacement','Angular velocity')
axis([0 20 -4 4])

%Animation for a pendulum

bob = results(:,1);
ct = 1; 
for i= 1:length(bob)
    thetax=bob(i);
    % Initial coordinates
    xo = 0;
    y0 = 0;

    % Final coordinates
    x1 = l*sin(thetax);
    y1 = -l*cos(thetax);

    %plotting
    figure(2)
    plot([-0.5 0.5],[0 0],'linewidth',3.5,'color','b');
    hold on
    plot([xo x1],[y0 y1],'linewidth',3.5,'color','r');
    hold on
    plot(x1,y1,'o','markers',20,'markerfacecolor','r')
    hold off
    axis([-1 1 -1.5 1])
    pause (0.003)
   file_text=sprintf('damper%d.png',ct)
   saveas(gca,file_text);
   ct = ct+1;
 end   

Output :

1.

2. Link for the animation video:

https://youtu.be/RRf1TTof_Ec

Errors:

1. Error faced while plotting the figure due to the wrong sequence of the lsode function.

CONCLUSION :

1. From this program we have solved the equation of motion of a simple pendulum successfully.
2. We have plot a graph of time vs (angular displacement and angular velocity).
3. We have created an animation of simple pendulum.