Second Order ODE PENDULUM
In order to solve the 2nd Order Ordinary Differential Equation, we have to first break it in two 1st Order ODE. The following ODE represents the motion of a simple Pendulum with damping:
`d^2theta)/dt^2 + b/m (d theta)/dt +g/l sintheta = 0 -------------(1)
where, b= damping coefficeint,
m= mass of the Pendulum ball in Kg,
g= acceleration due to gravity in m/s^2,
l= Length of the Pendulum in m.
We have to break this 2nd Order ODE in to two 1st order ODEs.
Let, theta = theta1
(d theta)/dt = (d theta1)/dt = theta 2 ---------- (2)
(d^2 theta)/dt^2 = (d^2 theta1)/dt^2 = d/dt(d theta1)/dt = (d theta2)/dt --------- (3)
for theta = theta 1,
(d^2 theta1)/dt^2 + b/m (d theta1)/dt +g/l sin theta1 = 0 ------(4)
from (2),(3) and (4),
(d theta2)/dt +b/m theta2 + g/l sin theta1 = 0.
so for solving the original 2nd order ODE, we have created two first order ODEs:
(d theta1)/dt = theta2 ;
(d theta2)/dt = -b/m theta2 -g/l (sin theta1).
Now we have to use the same concept and write a Program in Matlab to solve the 2nd Order ODE and produce an animation of a simple Pendulum.
It will be done using an in-built solver of ODE called \"ode45\".
Our Program will have two parts namely, 1) Main code 2) function code.
function code will be called in the main code for solving the ODE.
PROGRAM
PART 1 : MAIN CODE
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);
%solve ODE
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span, theta_0);
%plotting
plot(t,results(:,1),\'linewidth\',2)
axis([0 30 -5 5])
hold on
plot(t,results(:,2),\'linewidth\',2)
ylabel(\'plot\')
xlabel(\'time\')
grid on
%Animation
figure
ct=1;
for i= 1:length(t_span)
x_initial = 0;
y_initial = 0;
x_final = l*sin(results(i,1));
y_final = -l*cos(results(i,1));
plot([-1 1],[0 0],\'linewidth\',2,\'color\',\'b\')
hold on
plot([0 0],[0 -1.5],\'linewidth\',2,\'color\',\'g\')
hold on
plot([x_initial x_final],[y_initial y_final],\'linewidth\',3,\'color\',\'k\')
hold on
plot(x_final,y_final,\'Marker\',\'O\', \'Markersize\', 20,\'MarkerFaceColor\',\'r\')
hold off
axis([-2 2 -2 1])
grid on
pause(0.005)
M(ct) = getframe(gcf);
ct = ct+1;
end
movie(M);
videofile = VideoWriter(\'pedulum_animattion_ode\',\'Uncompressed AVI\');
open(videofile);
writeVideo(videofile,M);
close(videofile);
PART 2: 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
this function code was called while solving ODE in the Main Code Program.
EXPLANATION
- First Command Window and Workspace are being cleared and closed.
- followed by defining the Inputs such as damping coefficient (b), acceleration due to gravity(g), mass of the pendulum ball(m) and length of the pendulum(l).
- Then Initial conditions were given theta[0;3], which are initial displacement ans initial velocity 3 rad/sec respectively.
- Time points is defined using linspace command from 0-20 seconds and taking 500 values between them.
- Now taking [t, results] as output array we solved for ODE using inbuilt Solver ode45 and we called the function ode_func as well and from thiswe will get angular displacement and angular velocity between the time span of 0-20 seconds.
- Using the Plot command we plotted the Angular displacement and angular velocity graph.
- Now using for loop we will generate the Animation of the Pendulum. For loop will run up to the total no. of values we are taking between the time span, in this perticular case 500 times.
- Now as the one end of the pendulum is fixed, so we start by taking x_initial and y_initial equals to 0 and calculate the final position of the pendulum as x_final and y_final.
- the cos component is taken as negative as the pendulum is moving in downward direction.
- now we plot a horizontal line (blue in color), which will act as fixed support.
- next we will plot a green vertical line which will align with the centre position of the pendulum.
- Now we will plot the string of the pendulum using the initial and final point cooridates i.e, x_initial , x_final and y_initial, y_final.
- next we will generate the pendulum ball using Marker command and setting the size at 20 and red in color.
- After this the Axis is defined followed by pausing the plot for 0.005 seconds.
- Next we get the individual frames using getframe command and save it in an Array called M.
- counter is incremented by 1
- using the individual frames a movie is made using movie command.
- a new AVI file is created where all the frames are stored using videoFile command.
- The same Videofile is opened, written and closed in the following steps.
Plot
In the above plot we can see the the Angular displacement (blue) and angular velocity(orange) with respect to the time span (0-20 sec).
ANIMATION
Errors in Programming
Error 1.
got this error because there was space between \"ode\" and \"45\" in line 18. resolved it by removing the space between them.
Error 2.
this error was resolved by adding a comma \",\" after \" \'linewidth\' \" in line 41.
Error 3.
Here you can see the pendulum ball isn\'t visible and the pendulum is oscilatting in upside down direction. this error was resolved by the taking the cos component as negative.
Error 4.
This error was resolved using \"hold off\" command so that we get the proper animation.