Solving second order ODEs

AIM: To write a MATLAB program to simulate the transient behaviour of a simple pendulum and to create an animation of its motion.

GOVERNING EQUATION:

The equation of motion of a simple pendulum with damping is-

`(d^2theta)/(dt^2) + b/m (d theta)/(dt) + g/Lsin(theta) = 0`

where,

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

L = length of the pendulum in `m`

m = mass of the ball in `kg`

b = damping coefficient in `(Ns)/m`

OBJECTIVES:

1. To simulate the transient behaviour at a simple pendulum.

2. To plot the graph of angular displacement `theta` and angular velocity `v` against time `t`.

3. To create an animation of its motion.

MATLAB CODE:

1. Main script

close all
clear all
clc
% inputs
b = 0.05;
g = 9.81;
l = 1;
m = 1;

% initial condition
theta_0 = [0; 3];

% time points
t_span = linspace(0, 20, 600);

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

% Plotting
figure(1)
plot(t, results(:, 1), 'linewidth', 2);
hold on;
plot(t, results(:, 2), 'linewidth', 2);
xlabel('Time', 'fontsize', 14);
ylabel('Plot', 'fontsize', 14);
legend('angular displacement', 'angular velocity');
grid on;

theta = results(:, 1);
x_pos = sin(theta);
y_pos = - cos(theta);
ct = 1;
for i = 1: length(theta)
    figure(2)
    plot([0 x_pos(i)], [0 y_pos(i)], 'linewidth', 3);
    hold on;
    plot(x_pos(i), y_pos(i), 'ko','MarkerFaceColor','k','Markersize', 20);
    plot([-2 2],[0 0], 'k-', 'linewidth', 5);
    axis([-2 2 -2 0.5]);
    title('Animation of a Simple Pendulum', 'fontsize', 16);
    grid on;
    pause(0.03)
    M(ct) = getframe(gcf);
    ct = ct +1;
    hold off;
end
%Animation Recording
movie(M)
videofile = VideoWriter('pendulum.avi');
open(videofile);
writeVideo(videofile,M);
close(videofile);

Given,

g = 9.81 `m/s^2`

L = 1 `m`

m = 1 `kg`

b = 0.05 `(Ns)/m`

`theta_ 0` = 0 radian

`v_0` = 3 `(rad)/s`

`t_0 = 0` sec

`t_(end) = 20` sec

1. linspace() is used to create an array of time having equal intervals in its elements.

2. ode45() solves the differential equation specified in ode_func file.

3. Plot() command is used to plot graphs.

4. xlabel() and ylabel() is used to label the axes.

5. axes() is used to give upper limit and lower limit of the axes, so that plot area does not change.

6. getframe(), captures the current axes as it appears on the screen as a movie frame.

7. movie(), plays the movie defined by a matrix whose columns are movie frames produced by getframe().

8. VideoWritter(), creates a video file.

9. open(), opens a specified file.

10. writeVideo(), writes video in the specified file.

11. close(), closes a specified file.

2. Function to solve the derivatives

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

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

end

Here, the 2nd order differential equation is converted into two 1st order differential equations.

`theta_1 = theta`

`theta_2 = (d theta_1)/dt = (d theta)/dt`

`(d theta_2)/dt = (d^2theta_1)/dt^2 = (d^2theta)/dt^2 = -b/m theta_2 - g/Lsin(theta_1)`

ERROR ENCOUNTERED:

Error 1:

This error has occurred because in the function file I have defined the function as

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

I have entered the extra part(, t_span, theta_0)) as I copied the function definition from the function calling line of the main script.

After removing the extra part(, t_span, theta_0)), this error is rectified.

Error 2:

This error has occurred because in axis function I have given the y value as 0.5 to -2 (axis([-2 2 0.5 -2])). But in MATLAB we have to first give the lower value, i.e -2. This error is rectified by giving y value as -2 to 0.5 (axis([-2 2 -2 0.5])).

ANIMATION:

Simple Pendulum Animation

GRAPH: