Week 3 - Solving second order ODEs

Aim: solving a second order ODE using matlab

Objective:

1. Solving a second order ODE that represents the equation of motion of a simple pendulum with damping
2. To animate the motion of the pendulum

Equation:

`(d^2theta)/(dt^2)+(b*d theta)/(m*dt)+(g*sintheta)/l=0`

Where,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

Here the given ODE is Second order ODE for the pendulum motion. The pendulum motion depends upon the Newton second law where it leads to forming the second order ODE.

Pendulum is an ideal model in which the material point of mass m is suspended on a weightless and inextensible string of length L. In this system, there are periodic oscillations, which can be regarded as a rotation of the pendulum about the axis O

Period of Oscillation of a Nonlinear Pendulum is described by the second order differential equation.

(d^2theta)/(dt^2)+(b*d theta)/(m*dt)+(g*sintheta)/l=0`

• The second order ODE is break down to first order ode and is solved.

Code:

clear all 
close all
clc

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

%initial condition as for pendulum(disp and velocity)
theta_0=[0;3]

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

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

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')
legend('displacement','velocity')
% Animation
ct = 1;
for i=1:length(t)
    
    %initial condition
    x0=0;
    y0=0;
    
    %final condition
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));
    
    figure(2)
    % creatingbeam
    plot([-1 1],[0 0],'linewidth',6,'color','r');
    axis([-1.5 1.5 -1.5 1.5]);
    hold on
    %creating line
    plot([x0 x1],[y0 y1],'linewidth',3,'color','b');
    
    %creating sphere
    plot(x1,y1,'o','markersize',20,'markerfacecolor','b');
    hold off
    pause(0.03)
    M(ct)=getframe(gcf);
    ct=ct+1;
end

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

function defination:

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

results and plot:

Animation vedio:

https://youtu.be/q2qa7NDEqPg

Explanation of code:

1.Inputs for pendulum are provided.

b=0.05 (damping coeff)

g=9.81 (gravity)

l=1 (length)

m=1 (mass)

2.initial conditions, displacement and velocity are provided as theta_0

theta_0=[0;3]

3. Integration is done w.r.t. time. Therefore time values are provided from 0-20 sec

t_span=linspace(0,20,500)

4. initial(x0,y0) and final coordinate(x1,y1) values are assigned for the pendulum

%initial condition

x0=0;

y0=0;

%final condition

x1=l*sin(results(i,1));

y1=-l*cos(results(i,1));

5. ode45 is an inbuild ODE function which is used to solve second order ODE and ode_func is a function call

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

6. ode_func is a function defination which gives the output of the second order ODE of a pendulum.This funtion breaks the second oder to first order ODE and provides the output as a column vector of displacement and velocity.

[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];

where,

dtheta1_dt= angular dispacement

dtheta2_dt= angular velocity

7. The first plot is created between time and (displacement, velocity)

8.animation is created by usinfg a for loop and length of for loop is the length of time

for i=1:length(t)

9. Fixed beam,line and sphere are created sing the plot command. Marker command is used in plot to create a sphere.

plot([-1 1],[0 0],'linewidth',6,'color','r'); % creatingbeam

plot([x0 x1],[y0 y1],'linewidth',3,'color','b'); %creating line

plot(x1,y1,'o','markersize',20,'markerfacecolor','b'); %creating sphere

10. A movie is created by noting the frame value of each theta valve w.r.t time using ‘getframe’ command and by using counter(ct). The value of counter gets incremented after each for loop.

ct = 1

ct=ct+1

11. a movie is created by using commands like ‘movie’,’ VideoWriter’,’ writeVideo’,’open’ and ‘close’ and the movie is saved as 'pendulum.avi'

Errors:

Error1: Inverted pendulum in figure2.

Solution : coerdinate y1 should be assigned as negative

y1=-l*cos(results(i,1));