Week 3 - Solving second order ODEs

Week 3 : Solving second oder ODE

The solution of the given probelm is attached below. Animation of the pendulum motion with respect to time.

https://drive.google.com/file/d/1esTwyOJGyxMCIgGX5JWJq6iniNv0WCCq/view?usp=sharing

In the given problem we are asked to find the oscillation in a simple pendulum. Osceillation in a simple pendulum indicates the transient system. The motion is dependent on space and time. These motions are represented using second order ordinary diffential equations. In the given problem we make use of a pendulum with a mass of 1 Kg, length 1m. Gravitational constant =9.81 m/s^2.

In the given problem we initialise the constants of the pendulum.

Gravitational contant =g

g=9.81;

Lenngth of the rod=l
l=1;

Mass of the marker attached to the pendulum=m
m=1;

We initalise the initial position and the initial velocity of the pendulum in an arrray names theta_0. First row inicate the initial position of 0 and the mean velocity 0f 3m/sec
theta_0=[0;3]

In the given problem we are asked to find the the position of the pendulum for 20 seconds. Hence we created an array named t_span. The array creates 500 values from 0 to 20seconds range.
t_span= linspace(0,20,500);

We make use of the ode45 function to convert second order ordinary differential equation into two first order euations.
[t, results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
ode45 is the function implements Runga-Kutta method with variable timestep of effecient computation. The function initIally reduces the given equation of motion into 2 first order equation.

The function ode45 redirects to another function named ode_func which is an ordinary function. The output of ode45 function is the time and the angle of each timestep. The function ode_func takes the value of time, theta, damping coefficient, gravitational constant, length of rod and mass of the ball as inputs

The division of second order equation into 2 first order equation is dine in the named function.

Thus we have find out two dtheta values which is used to find the angle.

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

We plot the graph of of the solved ordinary differential equation with respect to time
plot(t,results(:,1))
hold on
plot(t, results(:,2))

Name the xaxis and the yaxis of the drawn plot
ylabel('Plot')
xlabel('time')

Now we try o animate the motion of pendulum with time. For that we initially create a variable to make the count ct and assign the value 1
ct=1;

Createed a for loop of length of 500 timesteps which is the length of the row matrix of the variable results.
for i=1:length(results(:,1))

We initialise the fixed point of the pendulum about which the pendulum oscillates. At point x0,y0 the pendulum oscillates at the value of 0
x0=0;
y0=0;

At the timestep t the postion of pendulum varies along the x and y direction. The end point of the pendulum is assumed a value of x1,y1 which represent the position of pendulum for each time. The timestep t is given by the loop variale i
x1=-l*sin(results(i,1));
y1=-l*cos(results(i,1));

We create subplots of 4 and 5 divisions in the x and y coordinates

subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12,13,14]);

We create a bar of length 2 units length in the x direction.
plot([-1 1],[0 0],'linewidth',4,'color','k')
hold on

We plot the fixed end position of the pendulum at the position x0,y0 and the end point x1,y1 and join them with the color red
line([x0 x1],[y0 y1],'linewidth',2,'color','r')
hold on

We create a ball at the end of the pendulum (x1,y1) with the marker command of size 20
plot(x1,y1,'-o','markers',20,'markerfacecolor','r');
grid on;

We define the area of plot to be required.
axis ([-1.5 1.5 -2 0]);

hold off

A plot is created on with pendulum attached with a ball at the end. The plot created has to be saved as a PNG file. For that a variable calld Text is produced into which the plot of each time frame is saved as a fileneame of Pendulum1, pendulum 2...etc.

Text=sprintf("Pendulum%d.png",ct)

The filename pendulamfollowed by the plot number is given by the count ct as shown in the figure. The produced png is file is later saved in to the system using the saveas command.
saveas(1,Text)

After obtaining the required figue we increament the value of couunt number ct by1.
ct=ct+1

The loop continues to giv a series of png files which has to be converted into an gif file using imagemagick application. We open a termial and navigate to the source where all the png files are saved.

magick convert Pendulum*.png movi.gif

Since the processing of imagemagick application is time consuming we have reduced the time step by deleting excessive figures. We have obtained the required animation. The link to the animation is given below

https://drive.google.com/file/d/1esTwyOJGyxMCIgGX5JWJq6iniNv0WCCq/view?usp=sharing

Octave programme for find the solution of undamped pendulum

clear all
close all
clc

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

%initial condition
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')
close


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

subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12,13,14]);
plot([-1 1],[0 0],'linewidth',4,'color','k')
hold on
line([x0 x1],[y0 y1],'linewidth',2,'color','r')
hold on
plot(x1,y1,'-o','markers',20,'markerfacecolor','r');
grid on;
axis ([-1.5 1.5 -2  0]);
pause(0.0003)
hold off
  Text=sprintf("Pendulum%d.png",ct)
  saveas(1,Text)
  ct=ct+1;
  
  end