Week 3 - Solving second order ODEs

#ODE of the motion of a pendulam

#Inputs

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

Here, 'b' is the damping coefficient, 'g' is the tacceleration due to gravity, 'l' is the length of string by which ball is attached, 'm' is the mass of ball.

#Initial condition

theta_0 = [0;3];

This array saves the initial values of theta i.e. angle made by string with respect to the normal, and its derivative with respect to the time i.e. angular velocity in rad/sec.

#Time points

t_span = linspace(0,20,500);

An array, 't_span' saves the values of time from 0 to 20 seconds in 500 intervals.

#Solve ODE

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

To save the results of ODE of motion of a pendulam, [t, results] is used. In this, a function of choise among ODE solvers i.e. 'ode45' is used to solve derivatives of theta with respect to time. This integrates the system of differential equations y′=f(t,y) from t0 to tf with initial conditions y0. Here, t0=0, tf=20, y0=theta_0. Also, there is one more function i.e. ode_func which is an user defined function. The code for the same is given below,

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

Here, the inputs required to run this function, are t, theta, b, g, l, m. 'theta(1)' and 'theta(2)' are considered for simplifying ODE of the motion of a pendulam. 'dtheta_dt' is an array which stores the values of first derivatives of theta1 and theta2.

#Animation for the motion of a pendulam

ct=1;

Here, 'ct' is a counter which is used for compiling figures of the motion of a pendulam at an instant.

for i = 1:length(results(:,1))
    x0=0;
    y0=0;
    x1=1*sin(results(i,1));
    y1=-1*cos(results(i,1));

Here, 'for' loop is used. In this, loop counter variable 'i' takes instant values of first column of an array 'results'. This first column contains 'theta' values for the time period of 20 sec. To animate the motion, one need to know the location of ball at specified intervals from 0 to 20 sec. The location of the ball can be specified in x- and y co-ordinate.Here, fixed point is specified as (x0,y0). As string makes an angle of theta w.r.t. normal, one can convert polar co-ordinates i.e. 'l' and theta into cartesian co-ordinates by using trignometric functions sine and cosine. Same logic is used for x1.y1.

plot([-1 1],[0 0],'linewidth',4,'color','k');
hold on
line([x0 x1],[y0 y1],'linewidth',2,'color','b');
hold on
plot(x1,y1,'-o','markers',20,'markerfacecolor','r');

 The first 'plot' command, creates a line which represents a fixed surface to which string is attached. 'line' command connects x0,y0 to x1,y1. Second 'plot' command creates a sphere at x1,y1 which represents a ball attached to the string. 

axis([-1.5 1.5 -1.5 0.5])
title('Motion of a pendulam')

'axis' command fixes the x- and y co-ordinates of graph showing the motion of a pendulam. 'title' command gives title to the graph. 

    pause(0.0003);
    hold off
    P(ct)=getframe(gcf);
    ct=ct+1;
end

'pause' command determines the interval between consecutive figures. 'getframe' command stores the figures in P(ct). Also, the counter 'ct' keeps on adding till it reaches to the end of the loop. 

#Movie formation

movie(P)
videofile = VideoWriter('Pendulam in motion.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,P)
close(videofile)

'movie' command saves the animation in different video formats. This code does operations such as record, open, write and save.

The link for the animation is given below,

https://drive.google.com/drive/folders/1QNuhFzs7ccE2dJzDDW2reTJt_Awwz1iA?usp=sharing

#Errors

Some minor error were there such as in y1=-1*cos(results(i,1), I forgot to put '-' negative sign, so the motion of pendulam was inverted. Other than this, there was no any error which showed an error message in command window.