Week 3 - Solving second order ODEs

%second Order differential equation
close all
clear all
clc

% Variables 
% Damping coefficient 
b= 0.05
% length in meters
l= 1
% mass in kgs
m= 1
% gravity in m/s^2
g= 9.81

%initial condition
theta = [0;3]
%Time 
t_span= linspace(0,20,200)
% ODE equation
[t, results]= ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta);

%plotting
figure(1)
hold on
plot(t,results(:,1))
plot(t,results(:,2))
xlabel('time')
ylabel('angular displacement and angular velocity')
grid on

figure(2)
hold on
x0 = 0
y0 = 0
ct= 1

%for loop for animation of pendulum
for i= 1:length(t)
    x1= l*sin(results(i,1))
    y1= -l*cos(results(i,1))
    
    plot([x0,x1],[y0,y1],'color','m','linewidth',10)
    hold on
    plot([-1.5, 1.5],[0,0], 'linewidth',8)
    hold on
    plot(x1,y1,'O', 'markersize',40, 'markerfacecolor', 'b')
    
    axis([-1.5 1.5 -1.5 1.5])
    grid on
    hold off
    
    pause(0.05)
    M(ct)= getframe(gcf);
    ct= ct+1;
end

%animation
movie(M)
videofile= VideoWriter('pendulum.avi');
open(videofile)
writeVideo(videofile, M)
close(videofile)

As shown above , It is Matlab program to solve second order ODEs such as Pendulum.

• First we have defined all the variables which we are going to use in second order differential equation.
• Secondly, we have defined the initial condition of theta which is nothing but the initial value of angular displacement and angular velocity 0 and 3 respectively.
• Then we have defined the time array using linspace command which we will choose 200 values between 0 to 20s and accordingly we will get 200 values of angular displacement and angular velocity.
• After that we are using ode45 command to combine theta and time span to generate a plot for the pendulum.
• Here ode45 is a Matlab inbuilt function which is saying that ode is a function of (t,theta) and we have defined ode function ode_func(t,theta,b,g,l,m).
• Now ode_func(t,theta,b,g,l,m) is a different program of matlab which we create to solve the equation of pendulum.

  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 we have derivative of theta1 and theta2 for the pendulum.

• To create a plot, we are using (:) for the raws and 1 and 2 for columns because here result would be column vector, which is having 2 columns and 200 raws and we will get the values of dispacement in column 1 and values of angular velocity in column 2.
• To create an animation movie first we need to define the initial values of pendulum 0,0.
• Then we are defining end values of pendulum which are x1,y1. Here we have to use for loop in which we are using i for the raws which shows the length of t.
• Afterwards we are generating plots for length of pendulum, axis and size of ball accordingly.
• we are using Marker command to create a ball of the pendulum.
• At the end we are using movie command to generate an animation movie of the pendulum same as we used in the program of forward kinematics.

Animation link:https://drive.google.com/file/d/17TqR6rRkr5v3-sdCEEu32hRJpPcbGqpF/view?usp=sharing

Also i added some screenshots showing errors while programming codes.

1. To overcome the error i got which is in screenshot i replaced the positions of t_span and theta in the line 21.