Solving Second Order ODE in MATLAB

Aim: To write a program in Matlab that will simulate the pendulum motion.

Objective:

1. Solving the given second-order ODE to simulate a simple pendulum.

2. Creating a movie out of simple pendulum animation.

Description:

Second-order ODE model many physical problems from fluid mechanics, heat and mass transfer, rigid body dynamics and elasticity. We can better understand the real-world problems by solving-second order ODE. The solving the given ODE we can understand the equation of motion of a simple pendulum with damping. This has a lot of applications in vehicle dynamics where the suspension system can be modeled using a set of mass and we can solve the ODE and figure out how vehicle dynamics are going to be.

When we have a higher-order ODE we break it down in a set of the first-order ODE.

given second order ODE,

here,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b= damping coefficient

Solving the given equation-

if,

`theta=theta1`

then,

`(d theta)/(dt)=(d theta1)/(dt)`

`(d theta1)/(dt)=theta2`

`(d^2 theta)/(dt^2)=(d^2 theta1)/(dt^2)`

`=d/dt((d theta1)/(dt))`

`=d/dt(theta2)`

Given ODE,

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

or,

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

From eqn.1 & 2

`(d theta2)/(dt)+b/mtheta2+g/Lsin theta1=0`

We have created two set of first order ODE in order to solve the second-order ODE.

1st ODE,

`(d theta1)/(dt) = theta2`

2nd ODE,

`(d theta2)/(dt)= -b/mtheta2-g/Lsintheta1`

Writing down in matrix form-

`d/dt[[theta1],[theta2]]`

`[[theta2],[-b/mtheta2-g/Lsintheta1]]`

1.0 Function code for the ODE: (Based on the solved ODE)

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

theta_1 = theta(1)
theta_2 = theta(2)
dtheta1_dt = theta_2;
dtheta2_dt = -(b/m)*theta_2-(g/L)*sin(theta_1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end

1.1 Program:

%  Program to simulate the transient behaviour of a simple pendulum

clear all
close all
clc

b = 0.05;
g = 9.81;
L = 1;
m = 1;

%Initial condition

theta_0 = [0;3];

% Time points

t_span = linspace(0,20,500);

% Solve ODE

[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,L,m),t_span,theta_0);

%Plotting

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Angular Displacement and Angular Velocity')
xlabel('Time')
grid on
legend('Angular Displacement','Angluar Velocity')
title('Oscillation Graph of Simple Pendulum')

%Plotting Pendulum

figure
ct=1;
for i=1: length(results(:,1))
    
    x_fixed = 0;
    y_fixed = 0;
    x_freeend = 1*sin(results(i,1));
    y_freeend = -1*cos(results(i,1));
    figure(2)
    plot([-1 1],[0 0], 'linewidth',5,'color','r');
    axis([-1.5 1.5 -1.5 1.5]);
    line([x_fixed x_freeend],[y_fixed y_freeend],'linewidth',3,'color','b');
    hold on;
    plot(x_freeend,y_freeend,'o','markersize',20,'markerfacecolor','b','markeredgecolor','g');
    title('Transient behaviour of a simple pendulum')
    grid on
    pause (0.01);
    M(ct) = getframe(gcf)
    hold off
    ct=ct+1;
end

movie(M)
videofile=VideoWriter('pendulum.avi','uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Output:

Figure-

Animation that represents the position of the pendulum-

2. Explanation of Code:

Let us first understand the function code-

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

theta_1 = theta(1)
theta_2 = theta(2)
dtheta1_dt = theta_2;
dtheta2_dt = -(b/m)*theta_2-(g/L)*sin(theta_1);
dtheta_dt = [dtheta1_dt; dtheta2_dt];

end

A function is a code that can be re-used by giving the correct parameter in a different code.

• Ode_func is a separate MATLAB function which we have created. Inside that, we are defining theta 1 and theta 2 which is represented as theta. Theta is our solution array, and b, g, L, m are our parameters of a pendulum.
• theta1_dt is our first ODE, it is equal to theta 2. theta2_dt is our 2nd ODE which we have derived earlier.
• Now we are combining them as an array called dtheta_dt and separate them as a column matrix with ; operator.

Explaining main code-

%  Program to simulate the transient behaviour of a simple pendulum

clear all
close all
clc

b = 0.05;
g = 9.81;
L = 1;
m = 1;

%Initial condition

theta_0 = [0;3];

% Time points

t_span = linspace(0,20,500);

% Solve ODE

[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,L,m),t_span,theta_0);

Here,

• Clear all, close all, clc command used to delete the variables created in the workspace and command window before writing a new program.
• b,g,L,m are our given parameters of a simple pendulum.
• Initial condition (theta_0): Here, we are solving the ODE system. There is two first-order ODE. one is displacement and second is velocity. The initial displacement of the pendulum is 0 and the initial velocity is 3. Now we are integrating it to find the solution.
• t_span: To integrate we need to find the value of time. Here, t_span is our time array. It goes from 0 to 20 sec with 500 values within.
• t, results: Since we are solving two ODEs, the results will 2 column matrix. The first column refers to displacement and 2nd column refers to velocity. ode45 is an inbuilt solver in a MatLab. Its main inputs are t and theta. ODE is a function of t and theta. Then We are defining our function. Inside that, we have our parameters t, theta, b,g,L,m which I have explained earlier in function code. After that, we are providing the integration with the time array (t_span) which describes how far the integration needs to be done and then the initial condition starting with the time (theta_0).

%Plotting

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('Angular Displacement and Angular Velocity')
xlabel('Time')
grid on
legend('Angular Displacement','Angluar Velocity')
title('Oscillation Graph of Simple Pendulum')

• Using plot command to plot the graph with the given statement which is our result (t, results). Here we are plotting Angular displacements vs Angular velocity (:,1) and (:,2).
• Using hold on command to retain the properties of existing results incase of multiple output graphs.
• Using grid on command to enable grid in the output graph.
• xlabel and ylable commands to give name to the x and y axis.
• Using legend command to to figure the Angular displacement and Angular velocity line.
• Using title command to give title name to the graph.

%Plotting Pendulum

figure
ct=1;
for i=1: length(results(:,1))
    
    x_fixed = 0;
    y_fixed = 0;
    x_freeend = 1*sin(results(i,1));
    y_freeend = -1*cos(results(i,1));
    figure(2)
    plot([-1 1],[0 0], 'linewidth',5,'color','r');
    axis([-1.5 1.5 -1.5 1.5]);
    line([x_fixed x_freeend],[y_fixed y_freeend],'linewidth',3,'color','b');
    hold on;
    plot(x_freeend,y_freeend,'o','markersize',20,'markerfacecolor','b','markeredgecolor','g');
    title('Transient behaviour of a simple pendulum')
    grid on
    pause (0.01);
    M(ct) = getframe(gcf)
    hold off
    ct=ct+1;
end

• Using figure command to display the graphical output.
• Using "For" loop command to repeat the value of results which has 500 values and it stores only one value. and the value of i is going to change as per our ODE result which we have solved.. And each time value changes it executes a set of commands that are present in for loop which is ending with "end" command.
• x_fixed and y_fixed represent the fixed end of the pendulum.
• x_freeend and y_freeend represents the moving sphere which moves on the basis of ODE result.
• plot command used for plotting a straight horizontal line from -1 to 1.
• Showing axis parameters with axis command.
• Using line command we are creating a moving line. At the one end it is fixed at x and y area and at the other end it is moving as per solved ODE result with linewidth of 3 and color red.
• Creating a sphere using plot command with the size of 20 and color blue with green edge color.
• Giving title name with title command.
• Using pause command to set an interval at each position.
• Using an array M and adding the getframe(gcf) function to create (stitching) an animation of all the plots.

• Using Counter command (ct=ct+1) to get a single position of the sphere at a time and each and every time for loop is getting executed it is going to increment 1 time.

movie(M)
videofile=VideoWriter('pendulum.avi','uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

Converting the animation plot of the above code into a movie with .avi extension with the above code.

• At the end of for loop MATLAB is going to make a movie with all the frames with movie(M) command.
• Creating a video with the command videofile.
  
  
• Opening a video file with open(videofile) command and then write the captured frame with writevideo command and then closing the video by close(videofile) command.

3. Error Faced:

The moving line with attached sphere was showing on the top area of the line instead of down. Solved it by adding (-) sign on the y_freeend.