Week 3 - Solving second order ODEs
function code :
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 code :
% A PROGRAM TO SOLVE DIFFERENTIAL EQUATION
clear all % to clear workspace
close all % to close all the windows like plot
clc % to clear the commands.
% inputs
b = 0.1; % DAMPING COEFFICIENT
g = 9.81; % M/S^2
l = 1; % METER
m = 1; % KG
% initial values
theta_0 = [0;3]; % initial angular displacement & initial angular velocity
% time point
t_span = linspace(0,10,300); % time range of motion
% solve ode
[t, results] = ode45(@(t,theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);
% ode45 takes 3 inputs which are the function to be intigrate & time span &
% initial condotion
plot(t,results(:,1)) % plot time vs result of 1st column
hold on
plot(t,results(:,2)) %plot time vs 2nd column
title('Hormonic motion of pendulum')
legend('angular displacement','angular velocity')
ylabel('plot')
xlabel('time')
%for loop
ct = 1;
for i=1:length(results(:,1))
x_start = 0;
y_start = 0;
x_end = l*sin(results(i,1));
y_end = -l*cos(results(i,1));
% plotting the animation
figure (2)
plot([-1 1],[0 0],'linewidth',5) % ploting of base line
hold on
plot([x_start x_end],[y_start y_end],'linewidth',3) % ploting of string
hold on
plot (x_end, y_end,'o','markersize',25,'markerFacecolor','m') % ploting of bob
axis([-1.5 1.5 -1.5 1]) % points in between motion takes place
hold off
M(ct) = getframe(gcf);
ct=ct+1;
end
% animation of the pendulum
movie(M)
videofile=VideoWriter('pendulum_animation.avi','UNCOMPRESSED AVI');
open (videofile)
writeVideo(videofile,M)
close(videofile)
A REPORT ON DAMPED HORMONIC MOTION USING PENDULUM
AIM:-
To solve an ordinary differential equation (ODE) of second order using Matlab & plot an animation using information.
Theory:-
Pendulum is a mechanism where a bob is attached with string & string is mounted to a support. we pull the bob at one side and leaves then the pendulum starts performing a harmonic motion i.e. it starts moves to & fro in between two points. But due to weight of bob is motion velocity continuously decreases & at the end it stays still at bottom i.e. weight of bob is normal to gravitational force.
The above figure shows the simple pendulum & various forces acting on it.
The above equation is a 2nd order differential equation. We use this equation to get results of damped harmonic motion.
Main Program : -
% A PROGRAM TO SOLVE DIFFERENTIAL EQUATION
clear all % to clear workspace
close all % to close all the windows like plot
clc % to clear the commands.
% inputs
b = 0.1; % DAMPING COEFFICIENT
g = 9.81; % M/S^2
l = 1; % METER
m = 1; % KG
% initial values
theta_0 = [0;3]; % initial angular displacement & initial angular velocity
% time point
t_span = linspace(0,10,300); % time range of motion
% solve ode
[t, results] = ode45(@(t,theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);
% ode45 takes 3 inputs which are the function to be intigrate & time span &
% initial condotion
plot(t,results(:,1)) % plot time vs result of 1st column
hold on
plot(t,results(:,2)) %plot time vs 2nd column
title('Hormonic motion of pendulum')
legend('angular displacement','angular velocity')
ylabel('plot')
xlabel('time')
%for loop
ct = 1;
for i=1:length(results(:,1))
x_start = 0;
y_start = 0;
x_end = l*sin(results(i,1));
y_end = -l*cos(results(i,1));
% plotting the animation
figure (2)
plot([-1 1],[0 0],'linewidth',5) % ploting of base line
hold on
plot([x_start x_end],[y_start y_end],'linewidth',3) % ploting of string
hold on
plot (x_end, y_end,'o','markersize',25,'markerFacecolor','m') % ploting of bob
axis([-1.5 1.5 -1.5 1]) % points in between motion takes place
hold off
M(ct) = getframe(gcf);
ct=ct+1;
end
% animation of the pendulum
movie(M)
videofile=VideoWriter('pendulum_animation.avi','UNCOMPRESSED AVI');
open (videofile)
writeVideo(videofile,M)
close(videofile)
Function Program : -
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
DESCRIPTION :
Program explanation :
- Initial commands like,
- Clear all – removes data from the workspace
- Close all – deletes the current figure or specific figure.
- Clc – clear the command window.
- After this we collect the various inputs like damping coefficient (b), acceleration due to gravity (g), length of the spring (l) & mass of the body ( m).
- Next we are going to decide the initial condition of position & displacement i.e. angular displacement of sting (theta) & angular velocity. After this we decide the timespan for our motion.
- After this we use ode45 solver. This is a function in Mat-lab which solves partial or complete differential equation.
- The value for ode45 are, the function to be integrated, time span, initial condition.
- [t, results] –it means that the result that have been calculated are w.r.t. time.
- After this for plotting a graph we simply use plot command and on command coz we are getting two results which velocity & displacement.
- Next we assign a horizontal limit which shows that our rigid support [x0 x1]. Then we assign the values for bob position. After this we assign our axis [-1.5 1.5 -1.5 1] limit in between which the bob going to be oscillate.
- After this for animation we use get frame command and write video command.
Problem faced:
I faced an error for too long later that I came to know that those are case sensitive events hence I have take care of the alphabets that they are small or capital.
Conclusion:
The differential equation is solve using ode function & after that we plot the graph & our animation videos are shown here in below link.
Link : - https://youtu.be/H3P8OM4hx_I
Window showing my matlab screen.