Week 3 - Solving second order ODEs

Program code for plot of simple pendulum and its animation

clc;
close all;
clear all;
b=0.05;
g=9.81;
l=1;
m=1;
%initial condition for displacement and velocity
%displacement is 0 and intial velocity is 3 m/s
theta_0=[0;3];
%time points
t_span = linspace(0,20,200);
ct=1;
%solve ode
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);%this is an 
%inbuilt function for getting the value for theta,angular velocity at each time interval.
figure(1);
plot(results)
figure(2);
    o=[0,0];
    for i=1:length(t);
        x=l*sin(results(i,1));
        y=-l*cos(results(i,1));
        P=[x,y];
        plot([-1 1],[0 0],[0,x],[0,y],'marker','o','markersize',3,'LineWidth',3);
        ball=viscircles(P,0.1);%it will create circle at the bob
        axis([-1.5 1.5 -1.5 1]);
        xlabel('x');
        ylabel('y');
        pause(0.003)
        M(ct)=getframe(gcf);
        ct=i+1;
    end
    movie(M);
    videofile=VideoWriter('pendulum motion.avi');
    open(videofile);
    writeVideo(videofile,M);
    close(videofile);

Animation for simple pendulum motion with damping :

Explaination :

In the above program we have created the data for theta and angular velocity for each time interval.t is divided into equally into 200 parts.for each time interval we have created theta and angular velocity,so that we can make a plot between time and theta and time and angular velocity.

Also we have writen code for motion of the simple pendulum using the data we have got from our ODE function(Ode_func()).using the data of angle theta for each time period we got x and y coordinate, wher x=sin(theta)and y=-cos(theta). The point P(x,y) is the one end of the rod, where bob of pendulum is fixed and (0,0) is the hinged point of the rod. plot shows the motion of a simple pendulum , where x=sin(theta) and y=-cos(theta).

Then we have written code for creating a movie.we have created a videofile named as'pendulum_motin.AVI'.

In this way we have sucessfully completed our program.

ODE FUNCTION

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

Plot for angle(theta) and angular velocity vs time :

Error faced while writing the program :

During writing the program i faced difficulty in writing the cde for motion of the pendulum.i didn't get how to write the code of oscilating motion. i got wrong plots.after trying by myself and took the help of internet also i failed to understand.After that i asked skill lynk support member,they told me how to write the code.i got the concept for writing the code.i understood the concept from 2R robotic arm animation program and successfully complete the writing of program code.