Week 3 - Solving second order ODEs

Aim:

Write a program in Matlab represents the equation of motion of a simple pendulum with damping

Theory:

The ordinary differential equations(ODE),are essential for the solving several analytical mechanics when a problem is involved with dependent and independent variables and derivatives. ode which consist of one or more functions of a independent variable and derivative of those function.

In partial differntial equation consists of more than one indepedent variable.

The second oder equation for simple pendulum is

In given equation

b = Damping coefficient = 0.06

g = Gravity(m/s2) = 9.81

l = length of the pendulum(m)=1

m=Mass of ball(kg)=0.1

Solving 2nd oder equation:

theta= theta1

in diffrentiate dtheta_dt=dtheta1_dt=theta2

d^2theta_dt^2= d^2theta1_dt^2 = dtheta2_dt

d^2theta_dt^2 = d^2theta1_dt^2=dtheta2_dt

d^2theta_dt^2+(b/m)dtheta_dt +(g/l)*sin(theta) = 0

d_theta2_dt+(b/m)theta2+(g/l)sintheta1 = 0

dtheta1_dt= theta2

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

Code:

clear all
close all
clc

%input variables

b =0.06;
g =9.81;
l =1;
m =0.1;
 
%initial condition
theta_0 =[0;5];

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

% ODE solve
[t,results] = ode45(@(t,theta) ode_func(t, theta, b,g,l,m),t_span, theta_0);
 
% Plotting of angular velocity and displacement
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('Time')
legend('Angular Displacement','Angular Velocity')

displ = results(:,1);

%Code for Pendulum

ct=1;


for    i=1:length(t);
    
    %Initial position
    x0 = 0;
    y0 = 0;
    
    %Final position
    x1 = 1*sin(displ(i))
    y1 = -1*cos(displ(i))
    
    figure(2);
    
    %line
    
    plot([-1 1],[0 0],'linewidth',3,'color','r')
    
    hold on
    
    %string
    plot([x0 x1],[y0 y1],'linewidth',3,'color','r')
    
    %Ball
    plot(x1,y1,'marker','o','markersize',25,'markerfacecolor','g')
    
    %Axis in plot
    
    axis([-3 3 -3 3])
    pause(0.02);
    hold off
    
    M(ct)=getframe(gcf);
    ct=ct+1;
end


movie(M);
videofile = VideoWriter('pendulum_oscillation.avi')
open(videofile);
writeVideo(videofile,M);
close(videofile);

Step Wise Explanation:

%Input variables

we need to add the parameters to the program so.

b =0.06;
g =9.81;
l =1;
m =0.1;

These are the variables needed to calculate the solution for the simple pendulum.

So then we need to find initial condition.so remember that we are solving a ODE system,and there are two ODE.

*velocity and Displacement

So that when solving 1st ODE we need to provide initial condition for the displacement.

%initial condition

theta_0 =[0;5];

theta_0 is assigned with initial angular displacement and angular velocity

so here we know that the displacement of the pendulum is 0 and initial velocity is 5m/s

Then we are going to integrating the ODE system to obtain a solution for which we need to find value for time .

t_span = linspace(0,20,500);

In the time span is the time array .it goes from 0 to 20 seconds in between the range of 500 values will be taken .

Then we have to call our ODE function .

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

Ode45 is the inbuild solvering in matlab.

(@(t,theta) --defines ode is a function of (t,theta) then basically say how ode function is defines or partial differential equation using numerical method

To calculate θ1 and θ2 ode command used in Matlab .By ode command, dθ1/dt and dθ2/dt is solved to get angular displacement and angular velocity . Here result in array with 2 coloums and row equal to number of time values .

% Plotting of angular velocity and displacement

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('Time')
legend('Angular Displacement','Angular Velocity')

Above command is used to plot the graph of Angular velocity and angular displacement.

displ=results(:,1);

This command is used to show displacement result ,that all 1st coloumn v alues of displacement are shown 500 values of displacement are given in the result.

Code for pendulum

%Code for Pendulum

ct=1;


for    i=1:length(t);
    
    %Initial position
    x0 = 0;
    y0 = 0;
    
    %Final position
    x1 = 1*sin(displ(i))
    y1 = -1*cos(displ(i))
    
    figure(2);
    
    %line
    
    plot([-1 1],[0 0],'linewidth',3,'color','r')
    
    hold on
    
    %string
    plot([x0 x1],[y0 y1],'linewidth',3,'color','r')
    
    %Ball
    plot(x1,y1,'marker','o','markersize',25,'markerfacecolor','g')
    
    %Axis in plot
    
    axis([-3 3 -3 3])
    pause(0.02);
    hold off
    
    M(ct)=getframe(gcf);
    ct=ct+1;
end


movie(M);
videofile = VideoWriter('pendulum_oscillation.avi')
open(videofile);
writeVideo(videofile,M);
close(videofile);

    

The above command is used to plot the pendulum animation ,In this postion of pendulum is defined and the pendulum size and ball size is defined. ball size is defined as Marker

The values are calculated as loop with the help of the command and all the animation results are stored in the array defined by (M)

Code To Call The Function:

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/1)*sin(theta1);
         dtheta_dt =[dtheta1_dt; dtheta2_dt];
         
end

Step wise explanation

The function is defined with the name of ode_func in which all the parameters from main program is called up in this function .The parameters are t,theta,b,g,l,m.

The second oder differntial equation is solved in to 1st oder DE in oder to calculate the value of theta1 and theta 2.

2 terms are stored in dtheta1_dt and dtheta2_dt.this value is going to be stored inside dtheta_dt in which values are transfer to the main project.

Result:

Plot diagram of Angular Displacement and Angular Velocity:

Error Faced While Doing Programme:

The error i had face is that i forget to give the code for displacement result . after i finding the error i make to correct the values in to clear the code.

Conclusion:

Differential equation of simple pendulum is solved using Matlab program .so we can say that solution obtained by program is correct and Hence we proved that.