Week 3 - Solving second order ODEs

% ode_func.m
% run this first standalone tab
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

% main_program.m

clc
close all
clear all

% Q1 : write a program for solving ODE

% Inputs
%damping_coefficient
b=0.05

%gravity meter/sec^2
g=9.81

%length_of_pendulum meter
l=1

 % mass kg
 m=1

 % initial condition
theta_0 = [0;3];

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

% solve ODE

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

figure("Name","Damping Motion of Pendulum")
hold on
% Displacement curve
plot(t,results(:,1),"LineWidth",3,"Color",'r')
hold on
% Velocity curve
plot(t,results(:,2),"LineWidth",3,"Color",'b')
xlabel('time')
ylabel('plot')
legend('Displacement(rad)','Velocity(rad/s)')
hold off
saveas(1,'Time_plot.png')
  
% Animation

theta = results(:,1);

figure("Name","Damping Oscillation of Pendulum")

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

% Rigid support 

% starting point of rigid support 
x2 = -0.5;
y2 = 1.5;

% Ending point of rigid support
x3 = 0.5;
y3 = 1.5;

% Input length
length_of_pendulum = 2

%initial condition hand Coordinates
    x0 = 0 ;
    y0 = 1.5;

for i = 1:length(theta)
    Theta1 = theta(i);

    % Thread length from initial point
    x1= length_of_pendulum*sin(Theta1);
    y1 = -length_of_pendulum*cos(Theta1);

    
    % rigid  support points
    plot([x2 x3],[y2 y3],'linewidth',6);
    hold on

    % Thread point from initial point to (x1,y1) point
    plot([x0 x1],[y0 y1],'LineWidth',4,'Color','g');
    hold on

    % To point Bob at (x1,y1)
    plot(x1,y1,'marker','o','MarkerSize',10,'Color','r','MarkerFaceColor','r');
%     figure dimmension
    axis([-4 4 -4 6])
    hold off
    pause(0.1)
    frame = getframe(gcf);
    writeVideo(videofile,frame)
    
end    
close(videofile);