Week 3 - Solving second order ODEs

Objective

The aim of the project is to write a program in MATLAB to simulate the motion of pendulum with given knowns:

• Gravity,g = 9.81 m⁄s^2
• Length of pendulum,L = 1 m
• Mass of ball, m = 1 kg
• Damping Coefficient, b = 0.05 kg/s
• Total Duration, t = 20 s

Intial Condition at t=0 s,

• Initial Displacement = 0 rad
• Initial Angular Velocity =3 rad/s

ODE: A classic example – Pendulum

In engineering, Ordinary differential equations ‘ODE’ is used to describe the transient behaviour of any system. Following the ODE explaining the motion of a pendulum

where,

g = gravity in m/s²

L = length of the pendulum in m

m = mass of the ball in kg

b=damping coefficient

The motion of pendulum depends on the Newtons second law. The above shown ODE describes the position of pendulum w.r.t time.

The ODE stated above is the second order differential equation. It can be solved more easily if the order of the equation can be reduced.

Reduction of Order

θ=θ1

dθ/dt=dθ1/dt=θ2

dθ2/dt=(-b/m)*θ2-(g/l)*sin(θ1)

so two first order equations are formed :

dθ1/dt=θ2

dθ2/dt=(-b/m)*θ2-(g/l)*sin(θ1)

MATLAB Solver: ode45

In order to solve the above reduced ODE, we require the time span and initial conditions.

The reduced ODE can be solved in MATLAB using inbuilt ODE solver, ode45 whose syntax is shown below:

[t,y] = ode45(odefun,tspan,y0): here y is the result and y0 is the initial condition

MATLAB CODE:

Explanation of code:

1. t_span is an array containing the time values from 0 to 20 subdivided into 200 divisions

2. [t,results]=ode45(@(t,theta) ans_func(t,theta,b,g,l,m), t_span ,[0;3]) in this line

• t=time outputs
• results contain two columns with various ouput angular displacement and angular velocity
• ans_func is the function
• t_span is the time limits
• [0;3] is the intial conditions

3. Results are plotted with plot command with adequate x label and y label

4. In order to find the position of the ball of the pendulum the final x and y position is calculated using 200 values

x1(i)=l*sin(results(i,1)) % x position of end points of pendulum
y1(i)=-l*cos(results(i,1)) %y position of end ponts of pendulum

5. plot command is used to plot the line plot([x0 x1(i)],[y0 y1(i)],'linewidth',2,'color','b') which indicates the length of pendulum

6. plot(x1(i),y1(i),'marker','o','MarkerFaceColor','r')---> This is used to construct a ball at the end position of the pendulum

7. ct is the counter to repeat the frames

8. M(ct)=getframe(gcf)--> M array stores all the frames generated

9. movie(M)- creating movies taking all the frames

clear all
close all
clc



% Solving the ODE that solves the equation of motion of simple pendulum 

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

%the time interval
t_span=linspace(0,20,200)% dividing the time into various subdivisions

%solving ode
[t,results]=ode45(@(t,theta) ans_func(t,theta,b,g,l,m), t_span ,[0;3]);%ode45 is used to get the values of theta and velocity at diffent times with initial angular displacement of 0 and initial velocity of 3


plot(t,results(:,1));% plotting time with theta values
hold on
plot(t,results(:,2));; % plotting time with velocity

xlabel('time')
ylabel('results')
legend('time vs theta','time vs velocity')

ct=1;

for i=1:length(results(:,1))
    
    x0=0;
    y0=0;
    
    x1(i)=l*sin(results(i,1)); % x position of end points of the pendulum
    y1(i)=-l*cos(results(i,1)); % y position of end points of the pendulum
    
    figure(2)
    
      plot(x1(i),y1(i),'marker','o','MarkerFaceColor','r'); %plotting the shape of ball of pendulum
      hold on
   
    plot([x0 x1(i)],[y0 y1(i)],'linewidth',2,'color','b'); %plotting the position of pendulum 
     plot([-1 1],[0,0],'linewidth',2,'color','k');
    axis([-2 2 -2 2]);
  
  hold off
    
    pause(0.1)
    M(ct)=getframe(gcf);
    ct=ct+1;
    
    
    
end
movie(M) %creating movie of all frames
videofile=VideoWriter('ode_program.avi','Uncompressed avi')
open(videofile)%opening the video
writeVideo(videofile,M) %Saving video file
close(videofile) %Closing the video file

Function used:

• Here inputs taken are t,theta,b,g,l and mass whose values are mentioned in the previous code
• Dtheta1 and Dtheta2 are the variables which stores the two state variables
• dtheta is a 2*1 array containing the Dtheta1 and Dtheta2 values
• Calling this function onto the ode45 command solves the two first order equations
• Name of this function is ans_func

function dtheta=ans_func(t,theta,b,g,l,m)
theta1=theta(1);
theta2=theta(2);

Dtheta1=theta2;
Dtheta2=(-b/m)*theta2 -(g/l)*sin(theta1);

dtheta=[Dtheta1;Dtheta2]
end

FIGURE OBTAINED:

OUTPUT OBTAINED MENTIONED IN WORKSPACE :

Youtube link for pendulum animation:

https://youtu.be/cQ1W3Lk6i9E