Week 3 - Solving second order ODEs

AIM

To calculate and plot second order equation by taking pendulum as an example.

To create an animation video of simple pendulum

THEORY

By taking damping coefficient as b,gravity as g,length of pendulum as l and mass of bob as m.where theta is the angle with respect to movemet of rope and t is the time of motion. we get the second order equation as

d^2(theta)/dt^2 = b/m*d(theta)/dt + g/l*sin(theta) = 0

CALCULATIONS

By considering theta = theta1 so d(theta)/dt becomes d(theta1)/dt which is equal to theta2

d^2(theta)/dt^2 = d(theta2)/dt

by substituting the values and rearranging the equation we get d(theta2)/dt = -b/m(theta2)-g/l sin(theta1).

To create animation plotting of pendulum we take the initial points as x0 and yo. the next points where the rope ends as x1,y1. then by trignometry we get the equation x1 = l*sin(theta) and y1 = -l*cos(theta).

PROGRAM

clear all
close all
clc

b = .3;
g = 9.81;
m = 1;
l = 1;

% initial conditions
theta_0 = [0,3];

%time points
t_span = linspace(0,20,100)

%solving ODE

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

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot');
xlabel('time');



theta = results(:,1);
c=1

for i = 1:length(theta)
    Theta = theta(i);
    
    % initial points
    x0 = 0;
    y0 = 0;
    
    %movement of pendulum
    
    x1 = l*sin(Theta);
    y1 = -l*cos(Theta);
    
    figure(2)
    % plotting string
    
    plot([x0 x1],[y0 y1],'color','g');
    hold on
    
    % plotting pendulum bob
    plot([-.6,.6],[0,0]);
    plot(x1,y1,'o','Markersize',15,'MarkerFacecolor','r');
    hold off
    axis([-1 1 -1.5 1]);
    pause(.03)
    M(c) = getframe(gcf);
    c = c+1;
end

  movie(M)
videofile = VideoWriter('pendulum','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

CODE EXPLANATIONS

The program is created in such a way which doesnt use more lengthy lines for doing differenciation. function made the differenciation easier

In the beginning the values of b,g,l and m are defined

in the next step theta_0 = (0,3) is defined where 0 is the initial displacement and 3 is the initial velocity.

t_span is the time for oscillation which is defined as linspace(0,20,100) wher time is 20 seconds and divided to 100 divisions.

this step is for solving second order equation that is

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

where ode45 is a function in matlab which solves partial differenntial equation using numerical methods.

so ode45 is the pivotal index of the program which solves the partial differential equation of simple pendulum.

here ode_func is a function program which is used for solving the equation by splitting it in to

dtheta1_dt = theta 2 and dtheta2_dt = -(b/m)*theta2-(g/l)*sin(theta1). and storing the values to

dtheta_dt = [dtheta_1_dt;dtheta2_dt] in which the values are stored in coloumn matrix.here the function program ends

then this function is called in the program by ode 45 which solves the differential equation.

in the next step the first figure of the velocity and didplacement is plotted by using the command

plot(t,[results(:,1)results(;,2)]. here the graph is produced and label x and y axis.

the next steps are for producing the animation of simple pendulum

it begins by defining theta = results(:,1) where colon is used to create all the rows of first coloum in matrix.

we create a for loop in order to create movement of rope and defines i = 1:length(theta) where i varies from 1 to the length of theta and the theta's each vale is stored in to Theta.

then initial points x0 and y0 are defined

movement of pendulum is calculated by x1 = l*sin(Theta) and y1 = -l*cos (Theta)

to plot the figure plot[x0 x1],[yo y1] is commanded and the hold on is used for plotting the bob

plot(x,y,'o','markersize','15','color','r','MarkerFaceColor','r') where markersize is the command for used to create a circle which is defined as 'o' in the code and MarkerFacecolor for coloring.

axis[-1 1 -1.5 0] is used to plot the graph inside the axis we require

M(c) = getframe(gcf) is for creating an animated video where Movie(M) is to store it as a movie

videofile is used to to store the animation as video format and then writevideo to acces the video file.

ERRORS

The error which i faced was to save it as avi format and when the small letteres of avi was changed to capital letters the errors was solved.

RESULTS

In the output the first column is displacement and the second is velocity

the velocity started with 2.82 m/s and decreased where as displacement started from 0 to .52

here is the plot graph

The animation video of a simple pendulum was produced