MATLAB Code for Solving ODE of Simple Pendulum

Aim:

• To write a MATLAB program that solves ODE which represents the equation of motion of a simple pendulum with damping.

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

• To write the MATLAB program that will stimulate the pendulum motions.

Theory:

In mathematics, an ordinary differential equation (or ODE) is a relation that contains functions of only one independent variable, and one or more of its derivatives with respect to that variable.

A simple example is Newton's second law of motion, which leads to the differential equation.

`m (d^2 x)/(dt^2 )=f(x)`

for the motion of a particle of mass m. In general, the force f depends upon the position of the particle x, and thus the unknown variable x appears on both sides of the differential equation, as is indicated in the notation f(x).

Second-order constant coefficient differential equations can be used to model spring mass systems.

An examination of the forces on spring mass system results in differential equation of the form

`mddotx+bdotx+kx=f(t)`

Where m represents the mass ,b is coefficient of damping force ,k is the spring constant ,and f(t) represents any net external forces on the system.

Similarly:

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

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

In the above equation,

g = gravity in m/`s^2`,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Damping is restraining of vibratory motion ,such as mechanical oscillation ,noise ,and alternating electric currents, by dissipation of energy. Unless a child keeps pumping a swing, its motion dies down because of damping. Shock absorbers in automobiles and carpet pads are examples of damping devices.

A damping coefficient is a material property that indicates whether a material will bounce back or return energy to a system. For example, a basketball has a low damping co efficient(a good bounce back). If the bounce is caused by an unwanted vibration or shock, a high damping coefficient in the material will diminish the response. It will swallow the energy and reduce the undesired reaction.

Solving 2^nd order ODE using MATLAB:

A typical approach to solving higher-order ordinary differential equations is to convert them to systems of first-order differential equations, and then solve those systems.

ODE of damped pendulum is 2^nd order ordinary differential equations , so we have to convert it into first order ODE.

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

In the above equation,

g = gravity in m/s²,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

`θ=θ_1` (equation 1)

From equation 1:

`(dθ)/(dt)=(dθ_1)/(dt)=θ_2 (equation 2)`

`(d^2 θ_1)/(dt^2 )=(d)/(dt) ((dθ_1)/(dt))=(dθ_2)/(dt) (equation 3)`

From equation number 2,3:

`(dθ_1)/dt=θ_2 ( equation 4)`

`(d^2 θ_1)/(dt^2 )+b/m (dθ_1)/dt+g/L sin⁡θ=0`

`(dθ_2)/(dt)+b/m θ_2+g/L sin⁡θ_1 =0`

so

`(dθ_2)/dt=-b/m θ_2-g/L sin⁡θ_1 ( equation 5)`

From equation number 4 and 5:

Writing equation in matrix format

`(d)/(dt) [(θ_1),(θ_2 )]=[(θ_2),(-b/m θ_2-g/2 sin⁡θ_1 )]`

After conversion use the MATLAB solver ode45 to solve the system.

The solver ode45 implements the Runge-Kutta(4,5) method. Such method is suited for solving ordinary differential equations by predictions.

[t,y] = ode45(@(t,y)odefun(t,y), tspan,y0), where tspan = [t0 tf], integrates the system of differential equations y′=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t.

Objectives:

1. Solve ODE which represents equation motion of simple pendulum with damping with help of MATLAB code.
2. Understanding of ode45 in MATLAB
3. Plotting the graph of variation of angular displacements vs time , variation of angular velocity vs time.
4. Over coming of errors if any.
5. Creating the animation plot of for motion of pendulum between motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

Body of content:

Given initial condition:

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/`s^2`

L, m, b & g represents length of pendulum, mass of bob, damping co-efficient & gravity.

Initial conditions:

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

analysis time 0-20sec

Below code will solve the ODE , will give the values of angular displacement and angular velocity for the given time slot and will plot the graphs of variaton of angular displacement vs time and variation of angular speed vs time.

clear all
close all
clc
b=0.05;
g=9.81;
l=1;
m=0.1;
%intial condition 
theta_0=[0;5];% initial angular displacement =0; intial angular velociy =5 rad/sec
%time point
t_span=linspace(0,20,500);
%slove ode
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
plot(t,results(:,1));
hold on;
plot(t,results(:,2));
ylabel('amplitude');
xlabel('time');
title({'Variation Of Angular Displacement vs Time','&','Variation of Angular Velocity vs Time'});
legend('angular displacement in rad','angular displacement in rad/s');
axis([0 20 -5 5]);
grid on;

values of angular displacement and angular velocity for given time slot:

so we can see that first values of results show the initial value of angular displacement and angular velocity.

Plot of angular displacement vs time and angular velocity vs time:

From above graph it can be observed that oscillation of pendulum are gradually damped as we can observe that angulalar displacement and angular velocity was gradually decreasing with time.

Now the below code will helps in stimulating animation of damped pendulum:

Error:

i got the above error because the in my code i have taken x1=l*sin((results(:,1))) and y1=-l*cos((results(:,1))) which is will lead to error while ploting plot([x0,x1],[y0,y1],'linewidth',2,'color','r'); because of dimensional mismatch

so we have to consider x1 =l*sin((results(i,1))) and y1=-l*cos((results(i,1)))

code after correction:

clear all
close all
clc
b=0.05;
g=9.81;
l=1;
m=0.1;
%intial condition 
theta_0=[0;5];% initial angular displacement =0; intial angular velociy =5 rad/sec
%time point
t_span=linspace(0,20,500);
%slove ode
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
figure (1)
plot(t,results(:,1));
hold on;
plot(t,results(:,2));
ylabel('amplitude');
xlabel('time');
title({'Variation Of Angular Displacement vs Time','&','Variation of Angular Velocity vs Time'});
legend('angular displacement in rad','angular displacement in rad/s');
axis([0 20 -5 5]);
grid on;
%below code will helps us to animate the motion of simple pendulum
figure(2)
ct=1;
%for loop
for i=1:length(results(:,1))
    x0=0;
    y0=0;
    plot([-1 1],[0 0],'linewidth',4,'color','k');
    hold on
    %position of bob
    x1=l*sin((results(i,1)));%(because (0,0) is fixed point ,x1 will be rightside of it during start)
    y1=-l*cos((results(i,1)));%(because (0,0) is fixed point ,y1 will be below of it during start)
    %plotting length of pendulum
    plot([x0,x1],[y0,y1],'linewidth',2,'color','r');
    %ploting of fixed end
    plot(x0,y0,'-s','MarkerEdgeColor','k','MarkerFaceColor','b','MarkerSize',10);
    %plotting of bob
    plot(x1,y1,'-o','MarkerEdgeColor','k','MarkerFaceColor','c','MarkerSize',20);
    hold off
    axis([-1.4 1.7 -1.5 0.5]);
    grid on
    titlestring=['Pendulum motion at time t=',num2str(t(i)),'seconds'];
    title(titlestring,'fontsize',10);
    
    M(ct)=getframe(gcf);
    ct=ct+1;
    
end
%code to animate the motion
movie(M)
videofile=VideoWriter('odesimplependulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

figure of animated plot:

Code regarding ode_func:

With help of below function we have converted second order eqation to first order:

function dtheta_dt = ode_func(t,theta,b,g,l,m)
theta1=theta(1);
theta2=theta(2);
dtheta1_dt=theta2;
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt];

Google drive link of animated plot:

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

Conclusion:

1. ODE has been solved using ode45 in MATLAB
2. Motion of simple pendulum has been animated.
3. It was observed that as damping co-efficient increses ,motion of pendulum gets damped faster
4. And from the above mentioned examples we can understand that damping has both advantages and disadvantages.