MATLAB program to solve second order ODE.

Objective :

To capture the motion of simple pendulum with damper by solving the governing equation for the motion which is basically a second order ODE and creating an animation using MATLAB.

The formula for 2nd Order ODE:

`(d^ 2 θ )/( d t ) ^2 + ( b/ m ) ⋅ (d θ) /(d t) + ( g/ L ) ⋅ sin ( θ ) = 0`

Code

1. First input parameters are assigned and then the initial conditions.
2. A time matrix is also assigned from 0 to 20s with 1000 points between them.
3. The time and results(angular position and velocity) are calculated using the ode45 function which takes input parameters such as the output of another function which calculates the differentiation of variable θ">θ w.r.t time and time matrix and initial conditions.

   function [dtheta_dt] = Pendulum_Damped_Code1(t,theta,b,m,g,l)
     
     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​

4. Angular position is plotted alongside with angular velocity.
5. Angular position vs angular velocity is plotted.
6. Using for loop and angular position the motion of the pendulum is animated.
7. The animation is saved then.

Error Faced and Solution :

1. 'results' variable is named as result in plot command. It has been then named correctly.

Result :

Graph :

https://drive.google.com/file/d/1jjBbr35An_LN6P1Ube0VhM7MTx9mA0DX/view?usp=drivesdk

Animation:

https://youtu.be/fRB8rubw0p8

Main Program:

clear all
close all
clc

%Input Variables

b = 0.05;            % damping coefficient
l = 1;               % length of pendulum(m)    
m = 1;               % Mass of bob(kg)
g = 9.81;            % gravity(m/s^2)

% Initial condition

theta_0 = [0;3];      % Angle is 0 and angular velocity is 5 m/s

% Integration limits

t_span = linspace(0,20,500);

% Solving ODE

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

% Plotting the results

figure
plot(t_span,results(:,1),'Color','r')
hold on
plot(t_span,results(:,2),'Color','b')
xlabel('Time')
ylabel('Position and Angular Velocity')
legend('Position','Angular Velocity')

% Plotting pendulum and capturing it as movie
figure
x_base = 0;
y_base = 0;
ct = 1;

for i = 1:length(t_span)
    x_end = l*sin(results(i,1));
    y_end = -l*cos(results(i,1));
    plot([x_base x_end],[y_base y_end],[-1 1],[0 0],'LineWidth',3)
    hold on
    plot(x_end,y_end,'Marker','o','MarkerSize',30,'MarkerFaceColor','r')
    axis([-1 1 -1.5 0.2])
    set(gca,'XTick',[], 'YTick', [])
    pause(0.03)
    hold off
    M(ct) = getframe(gcf);
    ct =ct+1;
end
movie(M);
videofile = VideoWriter('Simple_Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)