Week 3 - Solving second order ODEs

AIM:

1. Write the program that solves the ode of motion of simple pendulam with damping

2.Animate the motion of the pendulam

STEPS:

1. solving the Second order differential equation:

The given equation is in second order equation,to solve the it,the second order diff equation into first order diff equation.

above result is the ODE equation of simple pendulam.

2.PROGRAM FOR SOLVING SECOND ORDER DIFF EQN OF SIMPLE PENDULAM:

step1: program for ode function:

%ode function
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

EXPLANATION:

creating the function dtheta_dt and give the input arguments(t,theta,b,g,l,m) and write the solved equations as code and end the function.

step 2: program for solving second order ode of simple pendulam

clear all
close all
clc

%input
b = 0.05;
g=9.81;
l=1;
m=1;

%initial condition
ang_disp_vel = [0,3];

%time points
t_span = linspace(0,20,500);

%solving ode
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,ang_disp_vel);

figure(1);

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('angular displacement & angular velocity')
xlabel('time')
legend('angular displacement','angular velocity');

%creating animation of simple pendulam

ct=1;
for i = 1:length(results(:,1))
    ang_disp = results(i,1);
%initial condition
x0 = 0;
y0 = 0;

%displacement condition
x1 = l*sin(ang_disp);
y1 = -l*cos(ang_disp);

figure(2)

plot([-1 1],[0 0],'linewidth',2,'Color','r');
hold on
plot([x0 x1],[y0 y1],'LineWidth',3,'Color','b');
plot(x1 ,y1,'o','MarkerSize',15,'MarkerFaceColor','c');
axis ([-2 2 -2 2])
xlabel('x');
ylabel('y');
grid on
pause(0.003);
hold off
%movie creation
M(ct)= getframe(gcf);
ct = ct+1;
end

%movie creation 
movie(M);
videofile = VideoWriter('simple_pendulam','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile);

%END

EXPLANATION OF CODE:

STEP 1:

define the initial conditions for solving ode

length = 1 metre, damping coefficient = 0.05,mass of bob = 1 kg,gravity = 9.81

step 2:

define the initial conditions of angular displacement = 0 m and angular velocity = 3 m/s.

step 3:

give the time as time span using linspace command at equal intervals of (0 to 20 sec) .having 500 values.

step 4:

call the ode function to solve the ode equation for the given values.

[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,ang_disp_vel);

step 5: creating figure 1

create the figure 1 to plot the angular displacement and angular velocity changes with respect to time using plot command.

plot(t,results(:,1)) for angular displacement
plot(t,results(:,2)) for angular velocity

step 6: creating animation of simple pendulam

create the for loop and define the length of angular displacment results and define the string length of the pendulam

initial conditions x=0,y=0;

end conditions

x1 = l*sin(ang_disp);
y1 = -l*cos(ang_disp);

creating figure 2

figure 2 consist of animation of simple pendulam

plot the horizontal line,length of string and bob using plot commands

plot([-1 1],[0 0],'linewidth',2,'Color','r');----- horizontal line
hold on
plot([x0 x1],[y0 y1],'LineWidth',3,'Color','b');----length of string
plot(x1 ,y1,'o','MarkerSize',15,'MarkerFaceColor','c');------ bob size and color

step 7: creating the movie

movie is created by getting the current frames in movie[M] variables

M(ct)= getframe(gcf);

ct=ct+1 --------incremental value to get the images

and end the loop

step 8:

movie is write as a video file using videofile and writevideo command.

ANIMATION VIDEO LINK:

https://youtu.be/2SifS64D0vQ

ERRORS WHILE PROGRAMMING:

ERROR 1:

The pendulam was upside down forget to add the negative symbol in line 40

%displacement condition
x1 = l*sin(ang_disp);
y1 = -l*cos(ang_disp);---- line 40

Solution : Give negative symbol to change the direction.

ERROR 2:

PROBLEM: Animation of pendulam is not oscillated while creating a movie.

SOLUTION: Instead on sin and cos in end conditions I gave sind and cosd in degrees,then change it to rad by giving sin and cos only.

Result:

The program for solving second order equation of simple pendulam and animation of pendulam is created.

Thanks & Regards

B.Riyaz