PROGRAM FOR SOLVING SECOND ORDER ODE USING OCTAVE

AIM:-

• Write a program for solving second order ordinary differential equation.

INTRODUCTION:-

• In Engineering, ODE is used to describe the transient behaviour of a system. it helps to understand lot of important phenomenon about the system. A simple example is a pendulum.
• The way the pendulum moves depends on the Newton’s second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" with respect to time/
• Similarly, it can be used in vehicle dynamics where suspension system can be modelled as mass, spring and damper system. In this field by solving ODE we can get results about how mass is going to move with reference of time.

GOVERNING EQUATION:-

• We all know Newton’s laws of motion which we used to calculate the acceleration of an object, so that we could compute and draw its new position in the world. Now, we're going to use those laws to calculate the acceleration of a pendulum, to compute and draw its position along the swing.
• A pendulum is a bob suspended from a pivot. Obviously a real-world pendulum would live in a 3D space, but we’re going to look at a simpler scenario, a pendulum in a 2D space—the program canvas.
• In this case, however, the pendulum bob doesn’t simply fall to the ground, since it is attached by a string to the pivot point. Instead the bob swings towards its rest position, and we want to compute the acceleration of the bob along the swing. We'll call that its angular acceleration, since the pendulum is accelerating along an arc.

• This ODE represents the equation of motion of a simple pendulum with damping

`(d^2\theta)/(dt^2)+(b/m)*((d(\theta))/(dt))+(g/l)*(sin\theta)`..............(A`

When you have a higher order ODE to solve it you ned to break it down in to set of 1^st order ODE’s. So, Let's solve this eqaution

Take ,

`(\theta)=(\theta1)`

`(d(\theta))/(dt) = (d(\theta1))/(dt) = \theta2` ...............(1)`

`(d^2\theta)/(dt^2) = (d^2\theta1)/(dt^2) = d/dt((d(\theta1))/(dt)) = d/dt(\theta2)`..............(2`

so, let's eq(1) & eq(2) in eq(A)

`(d(\theta2))/(dt) + (b/m)*(\theta2)+(g/l)*sin\theta1`

`(d(\theta1))/(dt) = \theta2`.........(a)`

`(d(\theta2))/(dt) = -b/m*\theta2-g/l*sin\theta1`......(b)`

`d/dt[\theta1;\theta2] = [ \theta2;-b/m*\theta2-g/l*sin\theta1]` (combine both the equations)

where,

g = 9.81 (gravity in m/s2)

L = 1 (length of the pendulum in m)

m = 1 (mass of the ball in kg)

b=0.05 (damping coefficient)

OBJECTIVE:-

• Write a program that solves the above ODE.
• Program for generating animation of pendulum motion.
• Program for generating plot of velocity & displacement distribution with respect to time.

PROGRAM:-

• Main Program to achieve all objectives:-

% Program for soving second order ODE.

clear all
close all
clc

% input required for calculation
b = 0.05;                     % Damping co-efficient
g = 9.81;                     % Acceleration due to gravity (m/s^2)
l = 1 ;                       % Length of pendulum string (m)
m = 0.1;                      % Mass of Pendulum (kg)

% To solve ODE need to provide Initial condition for displacement
theta_0 = [0;5];

% For integrating the system to obtain solution need to define Time span
t_span = linspace(0,10,200);

% Call for ODE function
[result] = lsode('ode_func',theta_0,t_span);


% program for plot graph of velocity and displacement Vs time duration
figure(1);
hold on 
plot(t_span,result(:,1),'linewidth',3,'color','r')
plot(t_span,result(:,2),'linewidth',3,'color','g')
xlabel('time','fontsize',15);
ylabel('velocity/dispacement','fontsize',15);
legend('Angular displacement/tiime','Angular velocity/time')
title('velocity/dispacement distribution','fontsize',15)
grid on

% program for generate animation
result1 = result(:,1);
ct=1;
for i = 1:length(result1)
    RESULT = result1(i);
    
    x0 = 0;
    y0 = 0;
    x1 = l*sin(RESULT);
    y1 = -l*cos(RESULT);
    
    
    figure(2)
    % for reperesenting horizontal line
    plot([-1 1],[0 0],'linewidth',4,'color','g');
    hold on
    % for reperesenting point of connection between pendulum string to fixed point
    plot([0 0],[0 0],'o','markerfacecolor','b','markersize',10);
    hold on
    % for reperesent line between bob and fixed point(string)
    plot([x0 x1],[y0 y1],'linewidth',3,'color','b');
    hold on
    % for plotting bob at the end of string
    plot(x1,y1,'-o','markerfacecolor','r','markersize',30);
    hold off
    axis([-1.7 1.7 -1.7 1.7])
    grid on;
    pause(0.00003);
    title('Pendulum kinematics','fontsize',15)
    file_text =sprintf("pendulum_kinematics%d.png",ct)
    saveas (gca,file_text);
    ct = ct +1;
    
    end

PROGRAM Explaination:-

First of all we need to provide basic input parameters for calculation.

% Program for soving second order ODE.

clear all
close all
clc

% input required for calculation
b = 0.05;                     % Damping co-efficient
g = 9.81;                     % Acceleration due to gravity (m/s^2)
l = 1 ;                       % Length of pendulum string (m)
m = 0.1;                      % Mass of Pendulum (kg)

We have two first order ODE one id velocity (eq(a)) and second displacement (eq(b)). To solve these first order ODE we need to provide initial conitions.

% To solve ODE need to provide Initial condition for displacement
theta_0 = [0;5];

For integrating the ODE system to obtain a solution we need to define a time span.

% For integrating the system to obtain solution need to define Time span
t_span = linspace(0,10,200);

For calculation or solve ODE we need to provide equation which we derived, we will provide it by another program function just to make main program easy to understand and less messy.

function [dtheta_dt] = ode_func(theta,t,b,g,l,m)
  % input required for calculation
  b = 0.05;                     % Damping co-efficient
  g = 9.81;                     % Acceleration due to gravity (m/s^2)
  l = 1 ;                       % Length of pendulum (m)
  m = 0.1;                      % Mass of Pendulum (kg)
  
  theta1 = theta(1);
  theta2 = theta(2);
  
  dtheta1_dt =theta2;
  dtheta2_dt =-(b/m)*theta2-(g/l)*sin(theta1);
  dtheta_dt =[dtheta1_dt;dtheta2_dt];        % For combined them   
  end
  

Now, we have to call this function program in to our main program.we eill use lsode command. lsode needs the function name, the initial state vector, and the time to solve for.

% Call for ODE function
[result] = lsode('ode_func',theta_0,t_span);

from these data we can plot graph of velocity and displacement distribution with respect to time.

% program for plot graph of velocity and displacement Vs time duration
figure(1);
hold on 
plot(t_span,result(:,1),'linewidth',3,'color','r')
plot(t_span,result(:,2),'linewidth',3,'color','g')
xlabel('time','fontsize',15);
ylabel('velocity/dispacement','fontsize',15);
legend('Angular displacement/tiime','Angular velocity/time')
title('velocity/dispacement distribution','fontsize',15)
grid on

Result:-

in this by changing value value for damping co-efficient we can analys the behaviour of velocity and displacement distribution.

Program for generating animation of pendulum motion.

% program for generate animation
result1 = result(:,1);
ct=1;
for i = 1:length(result1)
    RESULT = result1(i);

in first line we name the variable result1 and we assign it all the values of result first column which is displacement, we will use for loop cammand and assign value of variable i ranging from 1 to length result1, which is desplacement values. then we assign result1(i) to variable RESULT.

    x0 = 0;
    y0 = 0;
    x1 = l*sin(RESULT);
    y1 = -l*cos(RESULT);

we provided the position of pendulum free an dfixed point.

Now commands for ploting graphs...

    figure(2)
    % for reperesenting horizontal line
    plot([-1 1],[0 0],'linewidth',4,'color','g');
    hold on
    % for reperesenting point of connection between pendulum string to fixed point
    plot([0 0],[0 0],'o','markerfacecolor','b','markersize',10);
    hold on
    % for reperesent line between bob and fixed point(string)
    plot([x0 x1],[y0 y1],'linewidth',3,'color','b');
    hold on
    % for plotting bob at the end of string
    plot(x1,y1,'-o','markerfacecolor','r','markersize',30);
    hold off
    axis([-1.7 1.7 -1.7 1.7])
    grid on;
    pause(0.00003);
    title('Pendulum kinematics','fontsize',15)

To get animation we plot graph for every value of displacement, save all the images and after stitching all of them together we will get animation. for this we will use 'saveas' command

to save the plot for every iteration, and 'sprintf' command for change the file name for every iteration. this will result in saving all the output files in the same directory.

Output:- Output animation video link :- https://youtu.be/fSSljgikI8o

ERRORS:-

1) i made typo, wrote colour inplace of color.

input:

-

Output:-

2) In saveas command i type (1, file_text), so i got fifty images of velocity distribition inplace of pendulum motion.which can solved by typing gca.

input:-

Output:-

CONCLUSION:-

• By increasing the value for damping co-efficient we can reduce the oscillation of the system.
• By this program we can solve ODE fro different systems and study th ebehaviour of the system.
• We can stydy the suspension system for vehicles.