Solving Second Order ODE by Using Matlab

Objective:

Write a Matlab Program to solve the second order ODE.

Given:

Here the given ODE is Second order ODE for the pendulum motion. The pendulum motion depends upon the Newton second law where it leads to forming the second order ODE.

Equation:

Conversion:

let theta = theta1 displacement

dtheta/dt=dtheta1/dt=theta2 angular velocity

dtheta2/dt = -b/m *theta2 + g/l *sin(theta1)=a

therefore , theta_dot =[ theta2;a]

So by numerically we can split this. In matlab everything we code as first order. By above it convert the second order ODE into two first order ODE.

To derive and calculate the value in matlab, a special command is used name as, ODE45.

syntax of ODE45:

[t,y] = ode45(@f(x),t,yo)

where t and yo are the initial condition.

Program:

clear all
close all
clc

% values
b=0.05; %damping co-efficient
g=9.81; %gravity
m=1; %mass
l=1; %length of the pendulum in metre.

%Initial condition
t_span=linspace(0,20,200); %time integrate
theta1=[0;3]; % angular dispacement =0, angular velocity = 3.

%solve differential equation
[t,results]= ode45(@(t,theta)pendulum_damped_ode(t,theta,l,g,b,m),t_span,theta1)

%plot the values
figure(1)
plot(t,results(:,1)) %time vs displacement
hold on
plot(t,results(:,2)) %time vs angular velocity
hold on
grid on
xlabel(\'time\')
ylabel(\'angular dispacement&angular velocity\')
legend(\'angular dispacement\',\'angular velocity\')

% Animation
ct = 1;
for i=1:length(results(:,1))
    
    %initial location
    x0=0;
    y0=0;
    
    %final location
    x1= l*sin(results(i,1));
    y1= -l*cos(results(i,1));
    
    figure(2)
    plot([-1 1],[0 0],\'linewidth\',6,\'color\',\'g\'); %creating a beam
    axis([-2 2 -2 2]);% fixed axis
    line([x0 x1],[y0 y1],\'linewidth\',2,\'color\',\'b\'); %creating the line from initial value to final value
    hold on;
    plot(x1,y1,\'o\',\'markersize\',20,\'markerfacecolor\',\'r\',\'markeredgecolor\',\'r\'); %create the sphere attached at end of the rope
    hold off;
    pause(0.1)
    
    %for anime
    z(ct) = getframe(gcf);
    ct = ct + 1;
end
%Animation converted
movie(z)
videofile = VideoWriter(\'pendulum.avi\',\'uncompressed AVI\');
open(videofile)
writeVideo(videofile,z)
close(videofile)

The above program is derived to calculate the second order ODE by using ODE45.

The main program contain the function calling which shows below,

Function Program,

function [dtheta_dt]= pendulum_damped_ode(t,theta,l,g,b,m)
theta1 = theta(1);%displacement
theta2 = theta(2);%velocity
dtheta1_dt = theta(2);%angular velocity
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);%acceleration
dtheta_dt =[dtheta1_dt;dtheta2_dt]
end
  %from ode45 it give result as 2 column, displacement and velocity because
  %it convert into numerical solution [t results]==== t spantime and
  %results is displacement and velocity respectively 
  

From the above program the ode45 solver calculate value for displacement and angular velocity with respect to the time.

From the above diagram it shows the size of the result where 2 column vector angular displacement and angular velocity respectively.

Graph:

Animation Video:

From the above graph and animation, it clearly shows that due to the low damping co efficient the motion is more and it is continuous. Increase the damping coefficient increase the resistance for the motion.

Requirement:

By knowing the syntax of the ODE45 and know the geometry and trignometric calculation we can easily code the program wihtout any error.

Conclusion:

From the above program, by using ODE45 we can simulate any higher order ODE which it convert into first order ode to run and get the results.