MATLAB function to solve an Ordinary Differential Equation (ODE) and simulate & animate the motion of a pendulum using the 'ode45' command.

I. AIM: To solve an Ordinary Differential Equation (ODE) and simulate & animate the motion of a pendulum using the 'ode45' command in MATLAB.

II. GOVERNING EQUATIONS:

`(d^2theta)/(d^2t) + b/m.(d theta)/dt + g/L.sin(theta) = 0`

where,

b = damping coefficient

m = mass of the ball in 'kg'

L = length of the pendulum (m)

g = acceleration due to gravity (`m/s^2`)

`theta` = angular displacement of pendulum

t = time

III. OBJECTIVES:

1. to understand how angular displacement changes with respect to time
2. to see how damping plays a role in the oscillations.

IV. PROGRAM

A.) A program to solve the above ODE and to plot the desired results

clear all
close all
clc

%% Inputs
m = 1;                                  % mass of the pendulum
b = 0.05;                               % damping constant
g = 9.81;                               % gravitational acceleration
l = 1;                                  % length of pendulum
t_span = linspace(0,20,500);            % time interval

%% Initial Conditions
theta_0 = [0; 3];                       % [Initial Angular displacement; Initial Angular velocity]

%% Solving the ODE
[time, results] = ode45(@(time, theta) ode_func(time, theta, b, g, l, m), t_span, theta_0);

%% Plotting
figure(1)
hold on
plot(time, results(:, 1), 'linewidth', 3)       % Plotting Angular displacement vs. time
hold on                                         
plot(time, results(:, 2), 'linewidth', 3)       % Plotting Angular velocity vs. time
xlabel('Time (s)')
ylabel('Angular Displacement (rad)')
legend('Displacement vs. time', 'Velocity vs. time')
hold on

%% Animation
hold on
figure(2)
ct = 1;
for i = 1:length(results)

    theta_p = results(i, 1);                    % converting the array of angles into a variable
    x0 = 0;                                     % Fixed end of the Pendulum
    y0 = 0;                                     % Fixed end of the Pendulum
    x1 = l*sin(theta_p);                       
    y1 = l*cos(theta_p);
    plot([x0 x1],[y0 y1], 'linewidth', 3, 'color', 'b');        % Plotting the various orientations of the pendulum
    line([-0.5, 0.5], [0 0], 'linewidth', 3, 'color', 'k')      % Plotting an arbitary fixed line
    hold on
    plot(x1, y1, 'o', 'markers', 20, 'markerfacecolor', 'm')    % Creating the ball
    hold off
    axis([-1 1 -0.1 1.5])
    grid on
    set(gca, 'Ydir', 'reverse')                                 % Reversing the Y-axis
    pause(0.01)
    M(ct) = getframe(gcf);
    ct = ct + 1;
    
end

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

1. We use MATLAB to conduct the following studies.
2. All the windows are closed and the workspace is cleared by using ‘close all’, ‘clear all’ and ‘clc’.
3. In the first section of the code we provide the inputs, such as mass of the ball, length of the pendulum, damping coefficient and the time interval.
4. Next we input the initial conditions and here the initial conditions are:
   • at time = 0, angular displacement = 0
   • and angular velocity = 3 rad/s
5. Now we solve the ODE function using the ode45 function.
6. The syntax of the ODE solver is [output1, output2] = ode45(@(independent variable, dependent variable) ODE_function, time span, initial conditions).
7. Here, we get two outputs, i.e, the time (s) and the desired results (angular displacement and angular velocity).
8. Hence, the output 'time' will be a column vector and the 'results' will be a vertical vector with two columns.
9. Now the angular displacement and angular velocity are plotted with respect to time.
10. Now to animate the pendulum we calculate all the [x, y] coordinates at every point of time. This is done by using the 'for' loop.
11. Now, we set the axis anf then we have to reverse the direction of the y-axis so as to maintain the proper orientation of the pendulum.
12. We load each frame in the variable 'M'.
13. Now we create a videofile using the 'VideoWriter' command and then open (using the 'open' command) it and write (using the 'writevideo' command) the vector M into the videofile and then we close it(using the 'close' command).

B.) A program to define the ODE function

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

theta1 = theta(1);
theta2 = theta(2);
dtheta_1 = theta2;
dtheta_2 = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta_1; dtheta_2];

end

1. We create a function 'ode_func' which takes in 6 inputs and returns one output.
2. The 6 inputs are :
   • the independent variable, time
   • the dependent variable, theta (angular displacement)
   • damping coefficient, b
   • acceleration due to gravity, g
   • length of pendulum, l
   • and mass of pendulum, m
3. Here, we convert the two degree homogenous equation into a one degree homogenous equation by substituting `(d theta_1)/dt = theta_2`.

V. PLOTS AND FIGURES:

Figure (1.): Angular displacement and angular velocity vs. time

VI. OUTPUTS:

The video can be found in this link:

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

VII. RESULTS:

1. We can see from figure (1.) that the angular displacement/velocity decreases with increase in time, this is because of the presence of damping.