Solving second order ODEs

Aim – Solving a second order Ordinary Differential Equation (ODE) using MATLAB.

Governing Equation –

The second order differential equation used for a simple pendulum is as following -

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

where,

b = damping coefficient = 0.05 Ns/m;

m = mass of bob = 1 kg;

g = gravitation acceleration = 9.81 m/s²;

L = length of the string = 1 m;

Above differential equation can be solved if we convert it into a first order differential equation.

Let,

`(dθ1)/dt = θ2`

and

`(d^2 θ1)/(dt^2 )= d/dt( (dθ1)/dt)= (dθ2)/dt`

then eq. becomes,

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

Now if we write it into matrix form then,

Objective – 1. To create a simulation for pendulum motion.

2. To write a program which can solve a second order Ordinary Differential Equation (ODE).

Code –

clear all
close all
clc

% Damping coefficient (Ns/m)
b = 0.05;

% Length of the string (m)
L = 1;

% Mass of bob (kg)
m = 1;

% Gravitational acceleration (m/s^2)
g = 9.81;

% Inputs ([theta at the start, velocity at the start])
theta_0 = [0;3];

% time points
t_span = linspace(0,20,1000);

%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))
xlabel('Time')
ylabel('Velocity(red) & Displacement(blue)')



% Simple Pendulum

ct = 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([x0 x1],[y0 y1],'linewidth',2,'color','b')
    hold on
    
    plot([-0.5 0.5],[0 0],'linewidth',4,'color','g')
    axis([-1.5 1.5 -1.5 1.5])
    hold on
    
    plot([x1],[y1],'o','MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','r')
    M(ct) = getframe(gcf);
    hold off
    ct = ct + 1;
end

movie(M)
v_f = VideoWriter('Simple_Pendulum.avi','Uncompressed AVI')
open(v_f)
writeVideo(v_f,M)
close(v_f)

Code for function ‘ode_func’ –

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

Explanation –

• First of all, we assign the input values of damping coefficient (b), gravitational acceleration (g), length of pendulum (L) and mass of the bob (m).
• Then we define a variable theta_0 which contains the initial angular displacement = 0 and initial angular velocity = 3. After that we assign the value of time span as 0 to 20 sec by the help of linspace command with 1000 intervals.
• Now to solve the ODE we use the command ‘ode45’ which is inbuilt and inside the command, we call the function ‘ode_func’ as calculations are done by the formula written in the function.
• Then we use the plot command to generate the plot between displacement and time and hold that ploy and superimpose another plot between velocity and time.
• Now to create an animation of pendulum, we have to use for loop with loop counter variable (i). Then we define the position of the end point of the string and plot it.
• Then to create the mass at one end of the pendulum, we use ‘markersize’, ‘markerfacecolor’ command and choose the size and color of the bob.
• Then to create a movie, we save all the frames in an array ‘M’ with the help of ‘getframe’ command and use hold off to create an animation. Then by increasing the counter (ct) by 1, we end the for loop.
• Now if we use movie(M), we will get the required simulation using command ‘VideoWriter’, we can dump all the frames into the v_f.
• Then we use ‘open’, ’writeVideo’ and ‘close’ command to run the program.

Result –

1. Displacement and Velocity plot with respect to time-

2. Screenshot of animation of Simple Pendulum –

3. Youtube link of animation –

https://youtu.be/rCyeVWF5XHw

Errors -

• Forgot to use the hold off command.
• After using the hold off command, error was rectified and required animation was obtained.

Conclusion –

Finally, the second order ODE is solved using MATLAB and the animation for a simple pendulum is created. By simply changing the inputs used, we can understand the effect of different parameters