Week 3 - Solving second order ODEs
Aim - To solve an ODE and use the calculated information to animate the motion of the pendulum
Introduction-:
In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum
The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.
Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.
Your objective is to write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.
Asimple pemdulum consits of a point mass suspended on a string oe wire that has negligible mass. When the bob is displaced through some angle pulling the string from its vetical position and is let go,it oscilates back and forth commencing the simple harmonic motion which is given by the following equation.
Calculation:
where-:
b= daming coeficient
m= mass of bob (in kg)
L= length of string (mt)
g= acceleration due to gravity(m/sec^2)
Following is the image of simple pendulum
Code explanation-:
Step 1: Providing the pre requisite values
FINAL STEP IS TO CREATE A MOVIE
clear all
close all
clc
%Giving the pre-requisite values for constants in equation
b = 0.05;
g = 9.81;
l = 1;
m = 0.1;
% Initial conditions
theta_0 = [0;5];
%time points
t_span = linspace(0,20,600);
%solve ode
[t, results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span,theta_0);
%plotting
figure(1)
hold on
plot(t, results(:,1))
plot(t, results(:,2))
title('Damping of pendulum harmonic motion')
legend('angular displacement','angular velocity')
ylabel('Plot')
xlabel('time')
%loop counter
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 animation
figure(2)
plot([-1 1],[0 0],'linewidth',4)
hold on
plot([x_start x_end],[y_start y_end],'linewidth',2)
hold on
plot(x_end,y_end,'o','MarkerSize',20,'MarkerFaceColor','auto')
title('Animation of damped oscilation of a simple pendulum')
axis([-1.5 1.5 -1.5 1])
pause(0.03)
hold off
M(ct) = getframe(gcf);
ct = ct+1;
end
%animation of the motion
movie(M)
videofile=VideoWriter('Pendulum Oscilation.avi', 'uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)
Program for function
function [dtheta_dt] = ode_func(~,theta,b,g,l,m)
theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
dtheta2_dt = -b/m*theta2 - (g/l)*sind(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];
end