Week 4.1 - Solving second order ODEs

Code and Result Explaination:

• Used ode45 to solve simple pendulum equation under damping.
• Function - 'ode_func' gives values of dtheta1/dt and dtheta2/dt into the ode 45 solver. this is done since the equation is 2nd order and needs to be reduced to 1st order ODE.
• Using for loop and we find the position of pendulum at each instance and create an animation for the same.
• Also, a plot is created between results (velocity and displacement) time which shows that velocity slowly decreases w.r.t time.

CODE:

1. Plot and animation code:

clear all
close all
clc

b = 0.05; %Damping Coefficient
m = 1; %mass
g = 9.81; %gravity
l = 1; %length 

%Initial condition for displacement and angular velocity
theta_0 = [0,3];
%time interval variable
timespan = linspace (0,20,1000);
%solve ode by using ode45
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),timespan,theta_0);

%plot for `
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Results (velocity and displacement)')
xlabel('Time')

%pendulum animation 
amp = results(:,1); % separate var for amplitude for using in animation
ct=1; %frame counter initialization
for i=1:length(amp)
    angle=amp(i);
    x0 = 0;
    y0 = 0;
    x1 = l*sin(angle); %position of pendulum w.r.t x axis
    y1 = -1*cos(angle); %position of pendulum w.r.t y axis

    %plot
    figure(2)
    plot([-0.5 0.5],[0 0],'linewidth',1.5); %celling line coordinate and plot 
    hold on %for animatiton we use hold on so it can be produced in the same plot
    plot([x0 x1],[y0 y1],'color','black','linewidth',1.5) % length (string of pendudlum) 
    plot(x1,y1,'o','markers',10,'markerfacecolor','red') %mass
    hold off 
    grid on 
    axis([-1.5 1.5 -1.5 0.5]) %axis definition
    pause(0.001)
    M(ct) = getframe(gcf); %frames of animation
    ct = ct+1;
end

%step 7: making movie
movie(M) %plays the movie frames
videofile = VideoWriter('simplependulum.avi','uncompressed avi'); %name and type of video
open(videofile)  %opens video file
writeVideo(videofile,M) %writes the movie m
close(videofile) %closes the video file
 

2. ODE function code:

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

    theta1 = theta(1)
    theta2 = theta(2)
    dtheta1_dt = theta2; % calculate dtheta1
    dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1); % calculate dtheta2
    dtheta_dt = [dtheta1_dt; dtheta2_dt];

end

PLOT:

Animation Link:

https://youtu.be/8ZdeslQ-kQY