Week 3 - Solving second order ODEs

AIM : To solve the second-order ode function using Matlab which depicts the motion of the pendulum and to create an animation of the same.

OBJECTIVE : To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping,

(d2θ/dt2)+b*m*(dθ/dt)+g*l*sinθ=d2θ/dt2+b*m*(dθ/dt)+g*l*sin⁡θ=0

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

GEOMETRIC PARAMETERS :

Length (L) = 1 metre,

Mass (m) = 1 kg,

Damping coiefficient (b) =0.05.

Acceleration due to gravity (g) =9.81 m/s2.

For Time period = 0-20 sec

Angular displacement=0, Angular velocity=3 rad/sec at time t=0

SUB CODE :

function dtheta_t = ode_func(t,theta,b,g,l,m)
theta1 = theta(1);
theta2 = theta(2);
dtheta1_t = theta2;
dtheta2_t = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_t = [dtheta1_t;dtheta2_t]
end

In this function program, the second-order ODE is simplified to the first-order ODE by introducing the term theta1 & theta2. Then theta1 & theta2 is differentiated and the dtheta_dt gives a column matrix.

MAIN CODE :

clear all
close all
clc
b=0.05;
m=1;
g=9.81;
l=1;
theta_0=[0;5]
t_span=linspace(0,20,200)
[t, results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m), t_span, theta_0)


figure(1)
plot(t,results(:,1),'linewidth',1)
hold on 
plot(t,results(:,2),'linewidth',1)
axis([0 25 -3 3])
xlabel('Time')
ylabel('Angular Velocity & Displacement')
legend('change in displacement','change in angular velocity')
hold off

figure(2)
plot(results(:,1),results(:,2))
xlabel('Position')
ylabel('Angular Velocity')


theta2 = results(:,1);
ct = 1;
for i = 1:length(theta2);
    theta3 = theta2(i);
    
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(theta3);
    y1 = -l*cos(theta3);
   
    figure(3)
    plot([-1.1 1.1],[0 0],'linewidth',2);
    hold on 
    plot([x0 x1],[y0 y1],'color','r','linewidth',2)
    hold on 
    plot(x1,y1,'o','markers',20,'markerfacecolor','b')
    hold off 
    grid on 
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.05)
end
plot(t,results(:,1),'linewidth',1)
hold on 
plot(t,results(:,2),'linewidth',1)
axis([0 25 -3 3])
xlabel('Time - sec')
ylabel('Angular Velocity & Displacement')
legend('change in displacement','change in angular velocity')
hold off

figure(2)
plot(results(:,1),results(:,2))
xlabel('Position - rad')
ylabel('Angular Velocity - rad/sec')


theta2 = results(:,1);
ct = 1;
for i = 1:length(theta2);
    theta3 = theta2(i);
    
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(theta3);
    y1 = -l*cos(theta3);
   
    figure(3)
    plot([-1.1 1.1],[0 0],'linewidth',2);
    hold on 
    plot([x0 x1],[y0 y1],'color','r','linewidth',2)
    hold on 
    plot(x1,y1,'o','markers',20,'markerfacecolor','b')
    hold off 
    grid on 
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.003)

    M(ct)=getframe(gcf)
    ct=ct+1;
end

movie(M)
videofile=VideoWriter('Pendulum', 'uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

In the above program, the fixed parameters are length, mass, damping coefficient, and acceleration due to gravity. Time and angle are given as the initial condition.In the tspan command, 200 values are taken in between 0 to 20 times period. Using the ode 45 command ODE function is introduced. Using the plot command velocity and displacement is plotted in the y-axis against time in the x-axis. Using the for loop command the mass at various points and various times are marked in the graph and are plotted. Using the movie command the plot is saved as an uncompressed avi file.

PLOTS :

ERRORS :

This error was solved by letting the movie complete for the total number of frames and not to close it while running.

MOVIE ANIMATION :

Link : https://drive.google.com/file/d/1V7gtH2x2CgS4NoV6P_tyB1CyOFUP0zNW/view?usp=sharing