Week 3 - Solving second order ODEs

Aim:

To write a MATLAB program to solve ODE which represents equation of motion of simple pendulum with damping.

Objectives:

• To Plot angular displacement and angular velocity as a function of time.
• To animate the motion of simple pendulum.

Theory:

Assumptions:

• Pendulum string is massless.
• Bob behaves like a point mass.
• The only resistance the system faces is the damping resistance.

Derivation:-

Here, equating the moments at hinge at the support-

`(d^2theta)/(dt^2)xxI = -(mgsintheta)l-(bv)l`

where,

I = M.I. of bob about hinge = `mxxl^2`

m = mass of bob (Kg)

l = length of string (m)

b = damping coefficient

v = tangential velocity of bob = `lomega`(rad/s)

dividing both sides by `mxxl^2` we get,

`(d^2theta)/(dt^2) + (b)/m((d theta)/dt) + g/lsintheta = 0`

for, `theta = theta1`

`(d^2theta1)/(dt^2) + (b)/m((d theta1)/dt) + g/lsintheta1 = 0`

and for, `(d theta1)/dt = theta2`

`(d theta2)/(dt) + (b)/m( theta2) + g/lsintheta1 = 0`

so, `[[(d theta1)/dt],[(d theta2)/dt]] = [ [theta2], [-b/mtheta2-g/lsintheta1] ]`

Body of Content:

To solve the ODE function, we need to define our function in separate program. As we can see we need input variables as b, g, l, m to be defined as well.

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); %dont use 'sind' here 'sin' only
dtheta_dt = [dtheta1_dt; dtheta2_dt];
end

Now Coming to main program, we already defined the input variables and the respective values. Initial conditions are specified as ang. displacement equal to 0 angular velocity equal to 3 and those values are stored in theta_0. Then to divide the time span into 400 steps we have used Linespace command. After that ode solver is directly implemented with the inputs from the function we defined earlier, this will give us angular displacement and angular velocity as two columns and rows being the time steps. And now we can use this data for plotting. in first plot we gonna plot angular displacement and angular velocity as a function of time. For this we need to use plot command twice first for angular displacement using first column and other for angular velocity using second column. For plotting2 we are actually animating the simple pendulum. First of all, we defined the hinge point using x0, y0. x1, y1 represents the other end where bob is attached. Using for loop we will get different values of x1, y1 for different timesteps, again plot command is used to trace complete trajectory of the bob. Using another plot command, we assigned the marker i.e. bob with ‘P’ corresponding to star shape and other variables like marker-size and marker face color. At the end of the for loop our standard format for creating a movie is used as shown in the program.

close all
clear all
clc

%Inputs
b = 0.05;   %damping coefficient
g = 9.81;
l = 1;
m = 0.1;

%initial condition
theta_0 = [0;3];   % displacement=0,angular velocity=3 at t=0

%time points
t_span = linspace(0, 20, 400); %dividing 20sec into 400 steps

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

%Plotting1
figure(1)
plot(t,results(:,1))    %First column of result until the end of rows
hold on                 %to retain data of first plot
plot(t,results(:,2))    %Second column of result until the end of rows
xlabel('time')
ylabel('plot')
legend('Angular Displacement','Angular Velocity')
grid on

%Plotting2
figure(2)
ct=1;
for i = 1:length(t_span)
    x0 = 0;
    y0 = 0;
    x1 = 1*sin(results(i,l));   %setting x-coordinates using angles in first column
    y1 = -1*cos(results(i,1));  %setting y-coordinates
    
    plot([x0 x1],[y0 y1],[-0.3 0.3],[0 0],'linewidth',2)    %for tracing the path and defining the support.
    hold on
    plot(x1, y1,'P','MarkerSize',25,'MarkerFacecolor','y')  %Assigning pendulum, 'P' corresponds to star-shaped marker.
    axis([-1.5 1.5 -1.5 1.5])   %Defining end limits af axes
    title('Animation of Simple Pendulum')
    grid on
    pause(0.00009)
    M(ct) = getframe(gcf);
    hold off
    ct = ct+1;
end
movie(M)
videofile = VideoWriter('Pendulum.avi','Uncompressed AVI'); %writing videofile name and it's type.
open(videofile)
writeVideo(videofile,M) %overwriting created videofile with the frames generated in 'for' loop.
close(videofile)

Final Output:

Plot1 is:

Here we can see that initial ang. displacement and ang. velocity at t=0 is what we already defined in the program. And it will further diminish with time because of damping.

Simulation:

Errors:

I forgot to include [] operator in axis command.