Week 3 - Solving second order ODEs

% ODE is used to describe the transient behavior of a pendulum

clear all
close all
clc

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

% initial condition
theta_0 = [0;3];

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

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

% plotting
subplot(5,4,[13,14,17,18]);
plot(t,results(:,1));
subplot(5,4,[15,16,19,20]);
plot(t,results(:,2));
plot(results(:,1),results(:,2));
xlabel('position radians');
ylabel('angular velocity');

ct = 1;
for i = 1:length(results(:,1))
    x = 0;
    y = 0;
    x1 = l*sin(results(i,1));
    y1 = -l*cos(results(i,1));
    subplot(5,4,[1:12]);
    plot([-1 1],[0 0],'linewidth',3,'color','b');
    hold on
    line([x x1],[y y1],'linewidth',3,'color','g');
    hold on
    plot(x1,y1,'o','markers',20,'markerfacecolor','r');
    grid on
    axis([-2 2 -1.5 0.2]);
    hold off
    G(ct) = getframe(gcf);
    ct = ct+1;
    
end
movie(G);
videofile = VideoWriter('pendulum.AVI','uncompressed avi');
open(videofile);
writeVideo(videofile,G);
close(videofile);

Function

% function used for ODE

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

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta(2);
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

Aim – Write a program in matlab that will simulate the pendulum motion.
Governing equation – d^2theta/dt^2 +b/m(dtheta/dt)+g/l(sintheta) = 0

Objective of the project – Objective is to write a program that solves the above ODE. This ODE represents the equation of motion of a simple pendulum with damping.
Body of the content - For simulation of a pendulum, we need input commands mentioned as below.         
 1) damping coefficient = b = 0.05
     Gravity = g = 9.81 m/s^2
     Length = l = 1m
     Mass = m = 1kg
     given as input. 
2) theta_0 = [0;3] is initial condition.
3) t_span = linspace(0,20,500) is 500 time points from 0 to 20 values .
4) Variable ‘dtheta_dt’ is defined for ode function contains inputs t, theta, b ,g, l, m.
5) theta1 = theta(1)
     theta2 = theta(2)
     dtheta1_dt = theta(2)
     dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1)
     2nd order differential eqn. converted to 1st order .
     dtheta_dt = [dtheta1_dt;dtheta2_dt]
     Solution vector/result stored in ‘dtheta_dt’.
     End is used for ending the function.
6) [t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
     Ode45 is inbuilt function .
     Ode is the function of t,theta.
     Ode function called with respect to the variable provided.
     The time for the integration to perform is t_span & theta_0 is initial condition where ‘0’ is 
     displacement & ‘5m/s’ is velocity.
7) subplot is used for plotting (5,4,[13,14,17,18])
8) plot is used for plotting time & results.
9) subplot is used for plotting (5,4,[15,16,19,20])
10) plot is used for plotting time & results.
11) plot(results(:,1),results(:,2))
       Plot is used for plotting of 1st & 2nd results.
12) xlabel is used to display position radians in x-axis .
       ylabel is used to display angular velocity in y-axis .
13) ct command is used for loop increment.
14) for loop used for multiple execution.
       x = 0 & y = 0 base 
       x1 = l*sin(results(i,1)) & y1 = -l*cos(results(i,1)) resulting position.
15) subplot is used for plotting (5,4,[1:12]).
16) Plot command is used for plotting x,x1 & y,y1 with line thickness 3 & blue colour.
17) Hold on is used for retaining the previous plot.

18) Line command is used for plotting line for pendulum with line thickness 3 & green colour.
19) Hold on is used for retaining the previous plot.
20) plot(x1,y1,'o','markers',20,'markerfacecolor','r') 
       Plot is used for creating pendulum. X1 & y1 is the end point of line where the pendulum to fix.
       ‘O’ is the shape of pendulum. Marker is the inbuilt command.
21) grid on is used for better display of dividing a page into major regions.
22) axis command is used to set x,x1 & y,y1 axis.
23) hold off sets the hold state to off so that new plots added to the axes clear existing plots and  
       reset all axes properties.
24) "G" is the variable in which all frames are merging for video.
25) after completion of one variable ct=ct+1 will increase.
26) end is use for ending for loop.
27) "G" variable is defined for movie in movie(G) . 
28) videofile command is used for naming the avi file & uncompressed format .
29) open command is used for opening the file.
30) writeVideo command is used for writing in "G" variable .
31) close is used for closing the video file .

Conclusion – Function is used as a separate program because of better usage. In function program we converted 2nd order differential equation to 1st  order .