Solving second order ODEs using MATLAB

Aim: Solving Second Order Differential equations using MATLAB.

Theory: A simple pendulum is one which can be considered to be a point mass suspended from a string or rod of negligible mass.

When displaced from its equilibrium point, the restoring force which brings it back to the centre is given by:

`F_"net"=−m⋅g⋅sinθ`

For small angles θ, we can use the approximation

`sinθ ≈ θ`

in which case Newton's 2nd law takes the form

`F_"net" = ma = −mgθ`

Even in this approximate case, the solution of the equation uses calculus and differential equations. The differential equation is

`(d^2 θ)/(dt^2 )+g/L θ=0`

However, the above case is assumed that there is no damping in the system. When we introduce damping, we get a modified equation as follows:

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

Where, b = damping factor and m = mass of the bob

To solve the above function, we use ode45 function in MATLAB. To know more about this function follow this link.

Objectives:

1. To write a program for second order ODE to show simple harmonic motion.
2. To show damping oscillation of simple pendulum in the form of graphical representation.
3. To animate motion of simple pendulum.

Program Outline

1. Clearing the previous all commands.
2. Setting up of inputs (b, g, L, m).
3. Setting initial conditions
4. Varying time span using linspace cmmand
5. Defining an ode function to store the above values of above matrix.
6. Using ode45 command to get the outputs
7. Using loop to get all the possible values of the problem
8. Animating all the figures to a video

Program Code

clear all
close all
clc

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

%INITIAL CONDITIONS
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);

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')
legend('displacement','velocity')

% IERATION USING FOR LOOP TO CREATE ANIMATION
ct=1;
theta_1=results(:,1);
for i=1:length(results(:,1))
    % INITIAL COORDINATES
    x_start=0;
    y_start=0;
    theta_2=theta_1(i);
    
    % FINAL COORDINATES
    x_end=l*sin(theta_2);
    y_end=-l*cos(theta_2);
    
    % PLOTTING FOR ANIMATION 
    figure(2)
    %Base of Pendulum
    plot([-1 1],[0 0],'linewidth',3,'color','black')                 
    hold on
    %String of Pendulum
    line([x_start x_end],[y_start y_end],'linewidth',5,'color','b')  
    hold on
    %Bob of Pendulum
    plot(x_end,y_end,'marker','O','markersize',20,'markerfacecolor','r')                     
    axis([-2 2 -2 1])
    title('Animation of Damped Vibrations in a simple pendulum')
    pause(0.000001)
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
end

% CREATING VIDEOFILE for M
movie(M);
% To create video file in avi format of animation
videofile=VideoWriter('Simple_Pendulum.avi','Uncompressed AVI');    
open(videofile);
writeVideo(videofile,M);
close(videofile);

Function Declaration

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];
end

Program Explaination:

1. Initially the workspace is cleared with all the unwanted variables.
2. Then inputs are defined for b, g, L and m.
3. Also the initial conditions along with the time pan is set as per the question stated.
4. Then a function is made names ode_func that takes six input arguments and returns an array of size 500 by 2.
5. Then a displacement v/s velocity plot is drawn using plot command
6. In order to aminate the pendulum, for loop is made which provides different locations of pendulums in order to make an animation.
7. Once the loop is executed, all these files are compiled into an animation using VideoWriter

Output Plot

Animation Video

Follow this link https://youtu.be/EuAQiMxM3zU

Errors Occured

• Improper function call

After encountering the error

Conclusion:

Second Order differential equations were developed using MATLAB. the necessary animations have also been displayed.

References

• http://hyperphysics.phy-astr.gsu.edu/hbase/pend.html#c4
• https://in.mathworks.com/help/matlab/ref/ode45.html?s_tid=srchtitle