Matlab program to solve a ODE representing the equation of motion of a simple pendulum with damping and simulating the motion of the pendulum.

The ordinary differential equation representing the equation of motion of a simple pendulum is

`(d^2theta)/dt^2+b/m*(d theta)/dt+g/lsintheta=0` --------------- eq 1

In the above equation,

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.

let us consider `theta=theta_1`

`implies(d theta)/dt=(d theta_1)/dt =theta_2`(say) ------------ eq 2

`implies(d^2 theta)/dt=(d^2 theta_1)/dt=(d theta_2)/dt`

the equation of motion can be writen as

`(d^2theta_1)/dt^2+b/m*(d theta_1)/dt+g/lsintheta_1=0`

`implies(d theta_2)/dt+b/m*theta_2+g/lsintheta_1=0`

`implies(d theta_2)/dt=-b/m*theta_2-g/lsintheta_1` -------------- eq 3

now we have two differential equations

`(d theta_1)/dt =theta_2`

`(d theta_2)/dt=-b/m*theta_2-g/lsintheta_1`

so in matrix from

`d/dt [[theta_1],[theta_2]] =[[theta_2],[-b/m theta_2-g/l sintheta_1]]` ---------------- eq 4

The MATLAB solver ode45 performs a direct numerical integration of a set of differential equations y' = f(t,y), y(t_o)=y_o from t_o to some final time t_f. The ode45 is a versatile ODE solver and we will use this solver in this project.

the syntax of ode45 is

[t,result] = ode45(@(t,theta) ode_fun(t,theta,l,m,b,g),t_span,theta_0)

theta_0 is an array of initial condition, ie. angular displacement and angular velocity at t=0,

t_span is an array of time points

ode_fun is the function defined as per equation 4

the solver gives us t and result

here the result is an array with two columns one with angular displacement and another with angular velocity.

now we can plot angular displacement vs time and angular velocity vs time

so the program will be

%Solving 2nd order ODE
clear all
close all
clc

%input
l = 1; % length of the string in m
m = 1; % mass of the bob in Kg
b = 0.05; %damping constant
g = 9.81; %acceleration due to gravity

%initaial condition
theta_0 = [0;3];

%time span
t_span = linspace(0,20,500);

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

%plot
figure(1)
plot(t,result(:,1))
title('Solution of ODE')
xlabel('time (s)')
ylabel('Angular displacement(rad)')

figure(2)
plot(t,result(:,2))
title('Solution of ODE')
xlabel('time (s)')
ylabel('Angular velocity(rad/s)')

ode_fun

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

Output -

now in order to animate the motion of the pendulum, we have to plot the position of the pendulum now to animate the string we will do simple mathematics, as we have already obtained the angular displacements for the time points, now as the one end of the pendulum is fixed so we are taking the co-ordinates of that end is (0,0) now for any instance of time let `phi` be the angular displacement so the coordinates of the another end of the string will be

`x = l*cos(phi-pi/2)` &

`y = l*sin(phi-pi/2)`

now why `(phi-pi/2)` because of the angular displacement obtained by solving the ode is measured along the y-axis and when we will plot it on Matlab the angle should be given with respect to the x-axis that is the reason for `(phi-pi/2)`.

Now in this way, we can plot the string but for the bob, we will use marker command.

Now we will run a loop from 1 to the length of the angular displacement array and for each angular displacement value we will plot the string and bob.

So as of now, we have generated the plots now we will use getframe function to capture the frames and store it in an array say M.

Now with the help of the movie command, we can generate the video animation of the motion of the spring.

so the Matlab code will be

%Solving 2nd order ODE
clear all
close all
clc

%input
l = 1; % length of the string in m
m = 1; % mass of the bob in Kg
b = 0.05; %damping constant
g = 9.81; %acceleration due to gravity

%initaial condition
theta_0 = [0;3];

%time span
t_span = linspace(0,20,500);

%solve ODE
[t,result] = ode45(@(t,theta) ode_fun(t,theta,l,m,b,g),t_span,theta_0);
ts=1;
%Animation
o = result(:,1);
for j = 1:length(o)
        phi = o(j);
        x = l*cos(phi-(pi/2));
        y = l*sin(phi-(pi/2));
        plot([0 x],[0 y],'-o','MarkerSize',20,'MarkerFaceColor','b')
        axis([-1 1 -1.5 0.5])
        title('Pendulum Animation')
        pause (0.03)
        M(ts)= getframe(gcf);
        ts=ts+1;
end
%video formation
movie(M)
videofile = VideoWriter('pendulum_animation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

ode_fun

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

output