Week 3 - Solving second order ODEs

Aim – Write a program to simulate pendulum motion.

Theory -

A simple pendulum consists of mass m hanging from a string length L and fixed at a pivot point. The differencial equation which represents the motion of a simple pendulum is-

`(d^2θ)/(dt^2) + b/m(dθ)/dt +g/L sinθ = 0`

m = Mass of the ball (Kg)

b = Damping coefficient

g = Gravitational acceleration (m/s^2)

L = Length of pendulum (m)

We are going to solve the second order ode by converting it into first order ode.

`(dθ_1)/dt = θ_2`

`(dθ_2)/dt = - b/m*θ_2 – g/l*sinθ_1`

input variables, we are going to use, are –

b = 0.05, L = 1 (m), m =1 (Kg), g = 9.81(m/s^2)

Code –

clear all
close all
clc

%input variables
b = 0.05;
g = 9.81;
l = 1;
m = 1;

%initial condition
theta_0 = [0;3];

%time points
t_s = linspace(0,20,300);

%solver ODE
[t,results] = ode45(@(t, theta) pendulum_damped_ode(t,theta,b,g,l,m),t_s, theta_0);

%plotting
figure(1)
plot(t, results(:,1))                % plot time vs angular displacement
hold on
plot(t, results(:,2))                % plot time vs angular velocity
xlabel("time");
ylabel("angular velocity and displacement");
ct=1;                                 % counter variable
for i = 1:length(results(:,1))
  disp = results(i,1);
  % coordinates of fixed end of string
  x0 = 0;                             
  y0 = 0;
  % coordinates of second end of string
  x1 = l*cos(disp-(pi/2));
  y1 = -l*sin(disp+(pi/2));
  % second figure window
  figure(2)
  plot([x0,x1],[y0,y1],'linewidth',3)                     % plot string
  hold on
  plot(x1,y1,'marker','.','markersize',40,'color','r')    % plot bob
  plot([-0.5,0.5],[0,0],'linewidth',8)                    % plot plateform where pendulum is attached
  hold off
  axis([-2,2,-2,2])
  pause(0.03)
  text = sprintf("pend%d.png",ct)                         
  saveas(2,text)                                          % saving each image
  ct = ct + 1;
end

Explanation of code –

1. First of all define the input variables under suitable name.
2. Define initial condition θ₀ = [0;3], a column vector.
3. Now use linspace command to define time span. Linspace command returns a row vector with linearly spaced elements.
4. Now apply ode45 function to solve ode. ode45 takes function handle, time interval over which ode will be evaluated and initial condition.
5. ode45 returns two output – a column vector which contains time and a matrix of unknowns.in this program it will return displacement and angular velocity as result.
6. Now plot the results we got from ode45 , time vs displacement and time vs angular velocity.
7. Set a counter variable and start a for loop which will run 1 to displacement values.
8. Set single value of displacement from its array, under a variable.
9. Set the coordinates of pivot point (first end of string) and string’s second end point.
10. Plot string, bob and a platform where string is attached with marker command.
11. Set axis and pause time.
12. Save each image of pendulum, changing through the loop to make a video
13. End the for loop.

Code of the function –

function [dtheta_dt] = pendulum_damped_ode(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];
endfunction

Conclusion -

we got a plot between angular displacement, angular velocity and time. intially angular velocity has value 3rad/sec. and angular displacement is starting with zero, as pendulum oscillates angular velocity and angular displacement decreases with time.

plot -

we can see pendulum animation using following link -

https://youtu.be/QDowvuY0N_w