Week 3 - Solving second order ODEs

Here We can simulate the transient behaviour of a simple pendulum and to create an animation of it\'s motion.

Below is the figure of simple pendulum :

Below is the Second order differential equation for simple pendulum:

Input Values:

b=0.05 (Damping coefficient)

g=9.81 (gravity in m/s^2)

l=1 (Length of pendulum in meter)

m=1 (Mass in Kg)

In intial condition Lets take,

Assume,

theta_0=[0;4];

Angular Dislplacement is 0 and Angular Velocity is 4

and

t_span=linspace(0,20,300);

Time span is 0 to 20 sec, with 300 intervals.

Now solve ODE,

In matlab ode45 is the inubuilt solver for ODE,

Here ODE is the function of t and theta , t and theta are dependent and independent variables.

Below is the syntax for ODE,

[t,result]=ode45(@(t,theta) ode_function(b,g,l,m,t,theta),t_span,theta_0)

We have called function of name ode_function into it and give a name to ode_function as [dtheta_dt]

First Simplify,

Assume,

theta=theta1................(1)

So, dtheta_dt=dtheta1_dt

dtheta1_dt=theta2.................. (2)

Then, d^2theta_dt^2=d^2theta1_dt^2

=d_dt(dtheta1_dt)

=d_dt(theta2)

d^2theta_dt^2=dtheta2_dt...........(3)

So, dtheta2_dt+(b/m)*theta2+(g/l)*sin(theta1)=0....(4)

dtheta1_dt=theta2........(5)

dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1).......(6)

Now we have two equations,lets club into matrix form.

d_dt[theta1; theta2]=[theta2; -(b/m)*theta2-(g/L)*sin(theta1)]

and dtheta_dt=[dtheta1_dt; dtheta2_dt]

To find dtheta_dt matrix, we will create ODE function

Note: theta is the solution array

Pelasefind below ODE function as name given \"ode_function\"

function [dtheta_dt]=ode_function(b,g,l,m,t,theta)
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

In above this function we have to find values for result. Result has two columns, 1st for displacement and 2nd velocity.

Thats why i am assigning,

theta1=theta(1) and theta2=theta(2)

The result part we will cover next,

Now, we will find values for displacement and velocity w.r.t. time in a matrix form of table result.

Then call above ode_function into main program,

close all
clear all
clc


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


theta_0=[0;4];

t_span=linspace(0,20,300);

[t,result]=ode45(@(t,theta) ode_function(b,g,l,m,t,theta),t_span,theta_0);

figure(1)
hold on
plot(t,result(:,1),\'color\',\'r\');
plot(t,result(:,2),\'color\',\'b\');

xlabel(\'Time\')
ylabel(\'Displacemt/Velocity\')

legend(\'Displacement\',\'Velocity\')

grid on

In ode45 synatx we defining value of b,g,l,m are assumed and then we define time span t_span and initial conditions by theta_0

The result of above progrm as belows:

I used hold on command to show the result into one graph or frame and legend command to differentiate displacement and velocity.

I got error when not use hold command;

In order to create animation vbideo i used counter (ct) command and movie command.

M is the variable to store all the output result figures and use for loop to complete time span length(t) i.e. 300 intervals.

close all
clear all
clc


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


theta_0=[0;4];

t_span=linspace(0,20,300);

[t,result]=ode45(@(t,theta) ode_function(b,g,l,m,t,theta),t_span,theta_0);

figure(1)
hold on
plot(t,result(:,1),\'color\',\'r\');
plot(t,result(:,2),\'color\',\'b\');

xlabel(\'Time\')
ylabel(\'Displacemt/Velocity\')

legend(\'Displacement\',\'Velocity\')

grid on


ct=1;

for i=1:length(t)
    
    
    x0=0;
    y0=0;
    
    x1=l*sin(result(i,1));
    y1=-l*cos(result(i,1));

    figure(2)
 
    plot([-2 2],[0 0],\'linewidth\',3,\'color\',\'r\');

    line([x0 x1],[y0 y1],\'linewidth\',4);
    
     axis([-2 2 -2 2])
     
     hold on
     plot(x1,y1,\'o\',\'markers\',20,\'markerfacecolor\',\'g\');
     hold off
     
     grid on
     pause(0.03);
     
     
     M(ct)=getframe(gcf)
     
     ct=ct+1;
end

movie(M)
videofile=VideoWriter(\'Second_order_ODE.avi\',\'uncompressed avi\')
open(videofile)
writeVideo(videofile,M)
close(videofile)

I got error when not fixing the reference axis :

in order to fix this i use axis([-2 2 -2 2])

I got another typographic error for wrongly used Videowriter command instead of used VideoWriter which.

Please find below link for Animation video of Simple Pendulum Single ODE

https://youtu.be/bPeRo_R44gM