Week 3 - Solving second order ODEs

Second order ordinary differential equation :-

An ordinary differential equation is a differential equation involving only one independent variable and combination of its derivatives. Such a differential equation whose order is 2 is called a second order ordinary differential equation.

The motion of a simple pendulum is governed by the ordinary differential equation of second order:

where

`theta`= angular displacement`

b=damping constant,

m=mass of ball in kg,

g=acceleration due to gravity in m/s2,

l=length of pendulum in m

As given in the problem statement, it is required to create an animation in MATLAB for the motion of a simple pendulum with the following data -

b=0.05,

g=9.8 m/s2,

l=1 m,

m= 1 kg

and initial conditions of angular displacement = 0 rad, angular velocity = 3 rad/s at time t=0 from t=0 to 20 seconds.

Here is the source code for the MATLAB animation :-

*********************************************************************************
MAIN FILE:-
*********************************************************************************
clear all
close all 
clc

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

%input initial conditions
theta_0=[0;3];

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

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

%plotting results
%title('Solution for ODE using ODE45')
%plot(t,results(:,1),'-o',t,results(:,2),'-o')
%hold ON;
x0=0;
y0=0;
disp1=results(:,2);
figure
%plotting for pendulum
for j=1:length(disp1)
    %clearing figure
    clf
    
    %instantaneous equations
    x1=l*cos(disp1(j)+ pi);
    y1=l*sin(disp1(j)+ pi);
    %plot instantaneous position
    plot([y0 y1],[x0 x1],'co',[y0 y1],[x0 x1],'b-','linewidth',2,'Markersize',10)
    axis([-1.2 1.2 -1.5 0.9])
    hold ON
    plot([-0.5,0.5],[0,0],'--r');
    
    %Animating for pendulum
    pendulum_ode(j)=getframe;
end

%save the movie
mywriter= VideoWriter('ODE_pendulum');
mywriter.FrameRate=10;

%opening,closing and writing movie
open(mywriter);
writeVideo(mywriter,pendulum_ode);
close(mywriter);.

*********************************************************************************
ODE Differential equation solver function handle:
*********************************************************************************

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

The code begins with

clear all

close all

clc

in order to clear all variables stored in the workspace, close all active plots and clear the command window.

The next part follows with the assignment of constants, namely, b,g,l and m, in the second order ODE.

Further, the code progresses on with

theta_0=[0;3];

t_span=linspace(0,20,500);

that gives the initial values of initial displacement of 0 radians and initial velocity of 3 rad/s and time stamp of observation of the simple harmonic motion of pendulum being studied from a time frame of 0 to 20 seconds with 500 discrete values.

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

ode45 is an ODE Solver in MATLAB and here, t and results are two matrices with t being a 500*1 matrix containing all time stamp values and results being a 500*2 matrix containing the corresponding values of angular displacement and angular velocities in the first and second rows respectively. Furthermore, @(t,theta) is the function handle being called for results with ode_func(t,theta,b,g,l,m) passing all arguments, t_span defining the time stamp or the integration interval and theta_0 gives the initial values.

ode_func is a function handle being called here with the definition

function [dtheta_dt]=ode_func(t,theta,b,g,l,m) defining function with arguments t,theta,b,g,l,m and dtheta_dt being the return value or matrix,

theta1= theta(1); assigning theta1 as the first column of theta matrix
theta2= theta(2); assigning theta1 as the first column of theta matrix
dtheta1_dt= theta2;
dtheta2_dt= -(b/m)*theta2-(g/l)*sin(theta1);
dtheta_dt= [dtheta1_dt; dtheta2_dt]; combining and converting into a 500*2 matrix of same sizes

For plotting the figure,

x0=0; defining points for point of pendulum pivot
y0=0;
disp1=results(:,2); fetching the angular velocity values
figure opens figure plot

For plotting pendulum,

for j=1:length(disp1) beginning of for loop

clf clearing figure

Instantaneous equations
x1=l*cos(disp1(j)+ pi); x component value of pendulum ball, pi is the suitable adjustment for time-interval
y1=l*sin(disp1(j)+ pi); y component value of pendulum ball

Plotting instantaneous position
plot([y0 y1],[x0 x1],'co',[y0 y1],[x0 x1],'b-','linewidth',2,'Markersize',10)

[y0 y1],[x0 x1] gives y and x values for plotting between dynamic value of ball and hinge point respectively, co : cyan colour and circles at plot points intended, blue solid line is intended to ne drawn and linewidth and markersize are chosen 2 and 10 respectively.

axis([-1.2 1.2 -1.5 0.9]) axes values have been suitably adjusted to get the best view
hold ON plots all the graphs on the subsequent time stamps on the same figure plot
plot([-0.5,0.5],[0,0],'--r'); horizontal red dashed line is intended to be drawn

pendulum_ode(j)=getframe; Animating for pendulum at jth instant by getting a frame and storing it in pendulum_ode matrix
end for loop ends here

%Saving the movie
mywriter= VideoWriter('ODE_pendulum'); creating an object and storing
mywriter.FrameRate=10; 10 frames per second has been chosen

%opening,closing and writing movie
open(mywriter);
writeVideo(mywriter,pendulum_ode);
close(mywriter);

***************************************************************

Errors faced while debugging :

There were numerous instances when final animation output didn't come out as desired. Out of them all, one instance is when the motion of the ball was observed to be as shown in the image crop.

Apparently, the gravitational force by earth is observed to be vertically upwards; the suitable modification that is rendered necessary is changing the time interval,i.e, increasing or decreasing the time interval by a value of pi radians.

%instantaneous equations BEFORE
x1=l*cos(disp1(j));

y1=l*sin(disp1(j));

%instantaneous equations AFTER
x1=l*cos(disp1(j)+pi);

y1=l*sin(disp1(j)+pi);

After changing time interval, the motion of pendulum is satisfactory.

Some error shown in command window is

Apparently, a "," was missing between the abscissa and ordinate points in line 36.

plot([y0 y1][x0 x1],'co',[y0 y1],[x0 x1],'b-','linewidth',2,'Markersize',10) was changed to

plot([y0 y1],[x0 x1],'co',[y0 y1],[x0 x1],'b-','linewidth',2,'Markersize',10) and the final results were desirable.

An AVI Video file has been uploaded on the link https://youtu.be/M41NoJx3I10 in YouTube.

A smoother and even more elegant animation has been uploaded on the link https://youtu.be/-MU3FqYGV7U in YouTube but with suitable modifications in the 'b' value(b=0.66).