Week 3 - Solving second order ODEs

clear all
close all
clc

%Inputs
g = 9.81; %gravity
l = 1; %lenght
m = 1; %mass
b = 0.05; %damping Coef 

%Initial Condition
theta_0 = [0;3];

%Time array
t_span = linspace(0,20,500);

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

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('plot') 


figure(2)
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([-1 1],[0 0],'linewidth',4,'color','k');
    hold on
    
    line([x0 x1],[y0 y1],'linewidth',2,'color','b');
    hold on; 
    plot(x1,y1,'-o','markers',20); 
    grid on;
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.003);
    hold off
    M(ct) = getframe(gcf)
    ct = ct + 1;   
end

movie(M)
videofile = VideoWriter('pendulum_motion.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

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

Objective:- Write a program in Matlab\Octave that will simulate the pendulum motion

Aim:- To simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0, to write a program that solves the following ODE `(d^2theta)/(dt^2)+b/m(d theta)/dt + g/l sin theta =0`

Theory:- 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 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 back and forth. This motion in the pendulum is given by the help of the ODE `(d^2theta)/(dt^2)+b/m(d theta)/dt + g/l sin theta =0`

Calculations :- For the pendulum we use the ODE thet defined as the equation of motion `(d^2theta)/(dt^2)+b/m(d theta)/dt + g/l sin theta =0`

for which the values where given as

g = gravity in m/s2 = 9.81 m/s2

L = length of the pendulum in m = 1 metre

m = mass of the ball in kg = 1 kg

b=damping coefficient = 0.05

For Simulation following things were provided;

Time Span = 0-20 sec,

Angular displacement = 0

Angular velocity=3 rad/sec

Time t=0.

Explaination:-

Step 1] To solve the first ODE the Inital conditions

theta_0 = [0;3];

this explains that the the initial displacement is 0 and the angular velocity is 3

Step 2] deifine the time array by using the linspace command

t_span = linspace(0,20,500);

as we are given the time as 20 sec. so linspace does it creates 500 point in between 0-20 or better yet divides 0-20 into 500 parts.

Step 3] Solving the ODE

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

ode45 is the ODE solver in Matlab

@(t,theta) we are saying or telling the these are the main inputs and they are the dependant and independant variables of the ODE basically oue ODE is a function of "t & theta".

ode_func(t,theta,b,g,l,m) this ia function created in matlab using the function command for getting t clean and neet to read code

t_span is the time array which tells us how far the integration will go on.

theta_0 these are the initial conditions

Step 4] Defining the function

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

here the theta is used as 1 in the solution array and hence in function it is defined as theta1 and theta2 for easier code to understand

dtheta1_dt = theta2; here we define `(d theta1) / dt = theta2` just for simplification in the problem

dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1); we just simplify the problem and get it in terms of dtheta2

dtheta_dt = [dtheta1_dt; dtheta2_dt]; here we just seperate it into a column matrix

Step 5] Plotting the graph

plot(t,results(:,1)) here we plot the first column which is displacement and we can know that as the value start from 0

hold on this command is given so the the graph of the displacemnet is held in the figure

plot(t,results(:,2)) here we plot the second column which is velocity and we can know that as the value start from 3

Step 6] Simulation

Simulation is done using the for loop so that every value can be taken into consideration one by one

so in the for loop

x0 = 0; here we define the stationery or the fixed end of the pendulum
y0 = 0;

x1 = l*sin(results(i,1)); here we define the moving end of pendulum using simple force vectors to calculate according to the theta
y1 = -l*cos(results(i,1)); here the '-' is present to get the pendulum 180 `degrees` so to get a proper simulation

plot([-1 1],[0 0],'linewidth',4,'color','k'); this is used to generate a line that will be acting as a wall to which the pendulum will be attached

line([x0 x1],[y0 y1],'linewidth',2,'color','b'); This is what creates the motion inside the simulation it is the line joining the stationery or the fixed end with moving end in the simulation

plot(x1,y1,'-o','markers',20); here we create a bob that is attached to the moving end and hece here the "x1 & y1" are used to plot so it moves with the string attached

M(ct) = getframe(gcf) now we create a frame for every step into the simultion for generation of the movie.

ct = ct + 1; this is the loop counter

Step 7] Creating a video file by using the videowriter in matlab.

videofile = VideoWriter('pendulum_motion.avi','Uncompressed AVI');

Results:-

Simulation:- https://drive.google.com/file/d/1EQJvrHo7Oe_7anGf8EDfMA8S9Gy2ldwF/view?usp=sharing