Week 3 - Solving second order ODEs

Aim:

To write matlab program to solve ode and to simulate pendulum motion.

Objective:

Here objective is to write a program to solve ode and to simulate the pendulum motion between 0-20second.

Introduction:

In the engineering ODE is used to describe the transient behaviour of a system. A simple pendulum is the best example. The way the pendulum moves depends on the Newton’s second law. When this law is written down we get a second order differential equation that describe the position of the Ball w.r.t. time. The differential equation is given as below:

`frac{d^2theta}{dt}+frac{b}{m}frac{d theta}{dt}+frac{g}{L}sintheta=0`

In the above equation g is the acceleration due to gravity, L is the length of pendulum, m is the mass of ball, b is the damping coefficient.

To solve ode it spilitted into 2 first order ode, the given ode is splitted as follows;

`frac{d theta}{dt}=frac{d theta1}{dt}=theta2`

`frac{d theta2}{dt}=frac{b}{m} theta2+frac{g}{L}sintheta1`

Explanation of code:

• In this program 1^st values of all the inputs are provided like l ,m, b, g. Then initial angular displacement and initial angular velocity are stored in an array.
• Then the time limit is set to 0-20 seconds with 500 values between them using linspace command.
• Then ode 45 function is defined as function of t and theta, after that ode function is defined which is a user defined function which takes t, theta, b, g, l, m as inputs, inside that we define theta1 and theta2, then we define 1^st order ode and 2^nd order ode, after that they are combined and put into an array, after that we provide time array and then inputs.
• The results are plotted from the solution of program and a plot is generated as output.

Explanation of animation of motion of pendulum code:

• In this program 1^st values of all the inputs are provided like l, m, b, g. Then initial angular displacement and initial angular velocity are stored in an array. Then the time limit is set to 0-20 seconds with 200 values between them using linspace command.
• Then ode 45 function is defined as function of t and theta, after that ode function is defined which is a user defined function which takes t, theta, b, g, l, m as inputs, inside that we define theta1 and theta2, then we define 1^st order ode and 2^nd order ode, after that they are combined and put into an array, after that we provide time array and then inputs.
• Then ode gives the results, where a array of 100 values is created in which 1^st column represents angular displacement(theta1) and 2^nd column angular velocity(theta2), after that theta is defined as the results of 1^st
• In order to create the animation initially counter is set to 1 and as the position of pendulum changes with r.t.time, it should be captured and should be captured, to do that for loop is used and loop is executed 100 times.
• X&Y coordinates of system are given, pendulum is hinged at(0,0), from hinged support a line is drawn and a bob is created at coordinates (x1,y1) using marker command of size 8, colour yellow using marker size and MarkerFacecolour command.
• Axis is bound from -2to2 in x direction, -2to2 in y direction, then counter is incremented and for loop is ended. Then a movie is created using movie command as avi file.

ode_func function code:

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

code for solving ode:

clear all
close all
clc

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

%initial condition
%initial angular displacemnt=0|initial angular velocity=5
theta_o = [0;5];

%time interval from 0 to 20 seconds
t_span = linspace(0,20,500);

% to solve ODE
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_o);

plot(t,results(:,1))
hold on 
plot(t,results(:,2))
xlabel('plot')
ylabel('time')
title('variation of angular displacement and angular velocity w.r.t.time')
legend('angular displacement','angular velocity')
grid on

plot of angular displacement and velocity:

code animation of motion of pendulum:

clear all
close all
clc

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

%initial condition
%initial angular displacemnt=0|initial angular velocity=5
theta_o = [0;3];

%time interval from 0 to 20 seconds
t_span = linspace(0,20,200);

% to solve ODE
[t,results] = ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_o);

theta = results(:,1);
ct = 1;
for i=1:length(results(:,1))
    theta_a = theta(i);
    
    %co ordinates of x&y
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(theta_a);
    y1 = -l*cos(theta_a);
    
    plot([-0.3 0.3],[0 0],'color','r','linewidth',6);
    text(0,0.2,'pivot');
    hold on
    
    plot([0 0],[-1.5 0.5],'k--','linewidth',0.5);
    plot([x0 x1],[y0 y1],'color','m');
    text(x1,y1+0.3,'string');
    hold on
    
    plot(x1,y1,'o','markersize',8,'color','g','MarkerFaceColor','g');
    text(x1-0.2,y1-0.1,'bob');
    hold off
    grid on
   
    axis([-2 2 -2 2]);
   
    M(ct)=getframe(gcf);
    ct = ct+1;
end

movie(M);
videofile = VideoWriter('pendulum animation.avi','uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Error:

Error is cannot find exact (case-sensitive) match for 'Movie'

Error here is due to wrong typing of command movie instead of typing movie, Movie has been typed, due to case sensitiveness of matlab it shown error. Then it has been resolved it by changing command as movie.

Conclusion:

In 1^st program a plot of 2 curves is created, where one curve represents angular displacement and another represents velocity.

In 2^nd program a animation of motion of pendulum is created.

Link for the animation:

https://youtu.be/9CA9G45go7Q