Week 3 - Solving second order ODEs

Objective: write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

Problem Statement:

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

use,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.

When you do this you will get the position of the pendulum as a function of time. Use this information to animate the motion of the pendulum.

• For each value in time
• Create an image that represents the position of the pendulum

Solution Procedure:

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of theses systems.

First of all convert second order ode into the two second order ode.

`d^2θ/dt^2+b/m.dθ/dt+g/lsinθ=0`

Where, b = damping coefficient

m = mass of the bob kg,

g = acceleration due to gravitym/s^2,

l = length of the string m.

`θ=θ1 `

`(dθ)/dt=(dθ1)/dt=θ2 -eq `

`d^2θ/dt^2=d^2(θ1)/dt^2=d/dt((dθ1)/dt)=d((θ2)/dt)−eq2 `

`(d^2θ)/dt^2+b/m.(dθ)/dt+g/lsinθ=0 `

`(d^2θ1)/dt^2+b/m.(dθ1)/dt+g/lsinθ1=0 `

`(dθ^2)/dt+b/mθ2+g/lsinθ1=0 `

`(dθ1)/dt=θ2 `

`(dθ2)/dt=−b/mθ2−g/lsinθ1 `

`d/dt[[θ1],[θ2]]=[[θ2], [(−bmθ2−glsinθ1)]] `

`d/dt[θ]=[[θ2], [(−bmθ2−glsinθ1)]`

Then assume theta as a solution vector.

So theta1=angular displacement

And theta2=angular velocity.

syntax for ode45:

[t,y] = ode45(odefun,tspan,y0),

where tspan = [t0 tf], integrates the system of differential equations y′=f(t,y) from t0 to tf with initial conditions y0. Each row in the solution array y corresponds to a value returned in column vector t.

matlab code for the solving second order ode is following:

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

clear all
close all
clc

%inputs 
b=0.05; %damping co.eff 
g = 9.81;  %gravity
l = 1;  %length of pendulum
m = 1;  %mass of ball

%initial conditions  angular displacement =0 and angular velocity=3
theta_0=[0;3];


%time points 0 to 20 sec
t_span=linspace(0,20,150);

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

%plot result 
figure(1)
plot(t, results(:,1),'g-')
hold on
plot(t,results(:,2),'r--')
hold on
xlabel('Plot')
ylabel('Time')
legend('angular displacement','angular velocity')

%Animation


counter=1;

for i=1:length(results(:,1))
    x0=0;
    y0=0;
    
    x1=l*sin(results(i,1));
    y1=-l*cos(results(i,1));
    
    figure(2)
    plot([-1 1],[0 0],'linewidth',10,'color','black')
    hold on
    plot([x0 x1],[y0 y1],'linewidth',1,'color','green')
    hold on
    plot(x1, y1,'o','MarkerSize',20,'MarkerFaceColor','r')
    
    axis([-1.5 1.5 -1.5 1.5])
    pause(0.0001)
    M(counter)=getframe(gcf)
    hold off
    counter=counter+1
   
end
%making the video file
movie(M)
videofile=VideoWriter('pendulum_ode_new.avi','Motion JPEG AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Code explanation :

First of all create a code which will solve the solution vector theta .(as explanied earlier)

then in second code give required input values.

use,

L=1 metre,

m=1 kg,

b=0.05.

g=9.81 m/s2.

select time span of 0 to 20 second using linspace command.

then use ode45 function which will give output as t and results.

then plot time vs.angular displacement and time vs. angular velocity in the same figure.

and for the animation take a counter variable and using required x and y coordinates plot the graph.

and for video file use getframe ,videowriter and writevideo function as shown in the code.

error while coding:

basic typo errors and during saving the videofile i have faced the error which is further solved by carefully go through the code.

Plot:

Animation :

youTube video link : https://youtu.be/4Jx3Aeso_jI

conclusion : The angular displacement and angular velocity are dying because of the damper.