Week 3 - Solving second order ODEs

Aim;

Aim: To write a program in MATLAB for describing the transient behavior of a system(A simple pendulum) by using 2nd order ODE’s.

Given data :

• g = 9.81 (gravity in m/s2)
• L =1 (length of the pendulum in m)
• m =1 (mass of the ball in kg)
• b=0.05 (damping coefficient)
• angular displacement=0
• angular velocity=3 rad/sec
• Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Solving the ODE'S

Given data:

`d^2 theta/dt^2+b/m d theta/dt + g/l sin theta`

In the above equation we are having the 2nd order ode.so, we have been convert into two first order ODE'S because the ODE45 accept only first order ode's.

theta = theta1

d theta/dt= d theta1/dt = theta2 ----->1

d^2 theta/dt ^2 = d^2theta1/dt = d/dt(d theta1/dt) = d/dt (theta2)

(from the give data) `^2 theta/dt^2+b/m d theta1/dt + g/l sin theta`

d^2 theta1/dt + b/m theta2 + g/l sin theta1 =0

d^2 theta1/dt = - b/m theta2 - g/l sin theta1

d theta/dt = - b/m theta2 - g/l sin theta1 ----->2

Arrange this two equation in the array form

`d/dt[[theta1],[theta2]] = [[theta2],[-(b/m) theta2 - (g/l) sin theta1]]`

Coding explanation

code-1

write a program to create a function then arrange this equation to find out the displacement and velocity with respect to time.

write a program and naming the function as ode_fun ensure the naming of function and saving the program as in folder name should be same (in line 2)

The value of theta(1) is to theta1 (in line 4)

The value of theta(2) is to theta2 (in line 7)

The value of dtheta1_dt is equal to the value of theta2 which we assumed in the above 1st equation

The value of dtheta2_dt is equal to the value of [`(b/m) theta2 - (g/l) sin theta1]`

creat an array of dtheta_dt of a 2-row and 1-coloumn matrix to show output of the equation.

end the function

so, in this program i can create a ode pendulam by calling this function.

% code 1
%creating the function which named as ode_func
function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

% The value of theta(1) is to theta1
theta1 =theta(1); 

% The value of theta(2) is to theta1
theta2 =theta(2); 

 % the value of dtheta1_dt is equal to the value of theta2
dtheta1_dt =theta2;  

 %the value of dtheta2_dt is equal to the value of -(b/m)*theta2-(g/l)*sin(theta1)
dtheta2_dt= -(b/m)*theta2-(g/l)*sin(theta1);

% creat an array of dtheta_dt of a 2-row and 1-coloumn matrix
dtheta_dt= [dtheta1_dt;dtheta2_dt];

end

code 2

In this program i can write a code to solve ODE function and and creation the animating video.

ODE function

The function of ODE 45 work with three inputs

1. Function has been integrated.
2. Time interval
3. Initial condition of displacement and velocity.

Once i provided the proper inputs means the ODE 45 should solve the equation.and its given the output for change in displacement and velocity of a pendulam with respect to time.

Loop for create animation

1. counter value ct=1 which mean to capture each and every position of the pendulam
2. using the for loop of i=1:t when the loop contineous in the ratio of time.
3. plot([x0,x1],[y0,y1],'linewidth',2) by plotting the co-ordinate of x0,x1 & y0,y1
4. Use the marker command to create a marker at the end of x1,y1 'marker size = 20'
5. stricting of all position of the figure in one array by using R(ct)=getframe(gcf) 'M' refers name of the array.
6. ct = ct+1 the counter of an image which adding of images getting from the loop one by one.

Creating the animation video:

1. Create a movie by using M(movie). "Mhttps://youtu.be/S8z205d71_Ahttps://youtu.be/-SoZImQV1lk" refers name of the array which having stiching of figures.
2. Create a empty video file and given name as 'ode pendulam' by using video writer command (videofile =VideoWriter('ode pendulam','Uncompressed AVI').
3. Open the video file by using open(videofile) command.
4. Then attach the stitching of all images to the empty video file by using writeVideo(videofile,M) command.
5. Close the video file by using the close(video) command.

code-2

% Solving second order ode pendulam

clear all
close all
clc

% Input data
b =0.05;     % Damping co-efficient
g =9.81;     % gravity in m/s2
l =1;        %length of the pendulum in m
m =1;        % mass of the ball in kg

% initial condition of displecement and velocity
% [:,1] represents as angular displacement & [:,2] represents as angular velocity
 theta_0 =[0;3];

%  time points
% time interval at 0 to 20 sec where 500 in between values are considered
t_interval =linspace(0,20,500);

% solve ode 
[t results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_interval,theta_0)

% plotting
figure(1)
hold on
% plotting a graph time vs result
angular_velocity=plot(results(:,2)); % plotting a graph of angular displacement
angular_displacement=plot(results(:,1)); % plotting a graph of angular velocity
xlabel('time'); % Named x axis as time
ylabel('angular-displacement & angular-velocity'); % Named y axis as angular velocity & angular displacement
legend('angular-displacement','angular-velocity');% highlight and differentiate the plotting x axis and y axis using legend command
% animation
ct=1 % assign the counter value as 1
% creating a loop for when the movement of pendulam with respect to time
for i= 1:length(t)
    x0=0; 
    y0=0;
    x1=l*sin(results(i));
    y1=-l*cos(results(i));
    figure(2)

    plot([x0,x1],[y0,y1],'linewidth',2) % plotting the co-ordinate of x0,x1 & y0,y1
    hold on
    plot(x1,y1,'o',MarkerSize=20); % use the marker command to create a marker at the end of x1,y1
    hold off
    axis([-3 3 -3 0]);
    pause(0.01)
    R(ct)=getframe(gcf); % stretch the all postion of figure's in one array by using getframe(gcf)
    ct=ct+1; % adding of counter plus one position
   
end
movie(R) % 'R' refers name of the array which having stiching of figures.
videofile =VideoWriter('ode pendulam','Uncompressed AVI');% create a empty video file and named as 'ode pendulam'.
open(videofile); % opening the created empty video file.
writeVideo(videofile,M); % adding the stretching of all images to make movie.
close(videofile) % close the video file.