Week 3 - Solving second order ODEs

%matlab code for 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


%matlab code 

clear all 
close all
clc

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);

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('vel & disp.')
xlabel('time')

theta = results(:,1);


%animation using for loop

ct = 1;
for i = 1:length(t_span)
    angle = theta(i);
    figure(2)
    title('simple pendulum')

    x0=0;
    y0=0;

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

    plot([x0,x1],[y0,y1],'linewidth',2)
    hold on
    plot(x1,y1,'.','markers',50)
    hold on
    plot([-3,3],[0,0],'lineWidth',3)
    hold off
    axis([-5 5 -2 2])
    pause(0.02)
    M(ct) = getframe(gcf);
    ct = ct + 1;
end


movie(M)
videofile = VideoWriter("simplependulum.avi","Uncompressed AVI");
open(videofile)
writeVideo(videofile,M)
close(videofile)

report:

->from the above code we can see that firstly the values of damping coefficient,gravitational acceleration,mass and length are given.

->then theta_0 = [0,3] which indicates that initial displacement and velocity is 0 and 3 respectively and a time span of 20 sec with 500 intervals.

->we use ODE45 which is an inbuilt matlab function to solve first order differential equations and we input our ode_func here.

->after solving the ode we get output as [t,results] where results are displacement and velocity and 't' is time and plot the result.

->for animation of pendulum we define the length and starting position of pendulum and to create the pendulum we use 'marker' command and give shape as '.' for it to be hollow.

->'for loop' is used to get the position of pendulum at each interval,length is changing at each angle and we run the for loop for theta=reesults(;1) for the entire first column and angle=theta (i) so we different displaacement value.

-> x1 and y1 to get the plots and we use pause command for time interval.

->to get all the frames we use getframe(gcf) and get motion of simple pendulum.

->to convert into a videofile we use 'movie' command we write the video and save it as avi file.

animation link:

https://www.youtube.com/watch?v=fpSBhPtd5SM