Solving second order ODEs in MATLAB

Aim:- To write a program that solves the following ODE which represents the equation of motion of a simple pendulum with damping. Additionally,Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

The ordinary differential equation is `(d^2 θ)/(dt^2 )+b/m*(dθ)/dt+g/l*sin⁡θ=0`

In the above equation,

g = gravity in meter per square second,

l = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

From the above equation, we create two ODEs because ode45 function can only solve first order ordinary differential equation.

`(dθ)/dt= θ_2,`

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

Given input values :-

L=1 metre,

m=1 kg,

b=0.05.

g=9.81`m⁄s^2`

Main Code:-

close all
clear all
clc

% Solving higher order ODE

% For a pendulum problem / damping system here.

b = 0.05;
m = 0.1;
l = 1;
g = 9.81;

% Initial condition
% look closely here its a column vector for as there are 2 ODE to solve
theta_0 = [0;3];

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

% solving ode
% here also the result vector has two set of solutions, so 2 plots below
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,m,l,g),t_span,theta_0);

plot(t,results(:,1),'color',[0 0.5 1])
hold on
plot(t,results(:,2),'color',[1 0.5 0])
xlabel('Time')
ylabel('Displacement(BLUE) or Velocity(RED)')
grid on
hold off

% Simulation
counter = 1;
for i = 1:length(results(:,1))
  x0 = 0;
  y0 = 0;
  x1 = l*sin(results(i,1));
  y1 = -l*cos(results(i,1));
  figure(2)
  plot([-0.5,0.5],[0,0], 'linewidth',3,'color','b');
  line([x0 x1], [y0 y1], 'linewidth',2,'color', 'r');
  hold on;
  
  plot(x1,y1, 'o', 'markers',18,'markerfacecolor','g');
  grid on
  axis([-1 1 -2 0.5])
  
  hold off
  M(counter) = getframe(gcf);
  
  counter = counter+1;
end

%creating Movie
movie(M)
videofile = VideoWriter('pendulum_animation.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function Code:-

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

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta2;
% Using degree(angle) here didnt yield my desired output check it yourself
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];
end

Explanation of the code:-

• The two ODE depends upon the time, theta and other parameters.
• In MATLAB “ode45” is an inbuilt solver, which is going to solve the ODE system.
• For animation film, we create an image for each value in the time that represents the position of the pendulum.
• For making bob the marker command is used.One end of the pendulum is fixed at (0,0) so that we can oscillate our pendulum.

Error Faced:

There was a syntax mistake here on line 57, Which was corrected by typing correctly "VideoWriter"

Outputs:

Here is a screenshot of the animation file,

The animation video link on youtube is this :- https://youtu.be/CU_d4qsSWXU

Conclusion :- In this way with the use of MATLAB study of transient motion of a simple pendulum with damping can be studied with the help of plots and animation, Similary we can solve other higher differnitial equations and its solutions & plots can be obtained easily.