Week 3 - Solving second order ODEs

simulation of pendulum

main code-

%g = gravity in m/s2,L = length of the pendulum in m,m = mass of the ball in kg,
%b damping coefficient.
b=0.2;
g=9.81;
l=1;
m=0.1;

% initial condition displacement=1,angular velocity=3rad/sec
theta_0=[0;3];

%time points
t_span=linspace(0,20,600);

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


ct=1;
%figure 1 using the movie command to create an animation of pendulum
figure(1)
for i=1:length(t_span)
    %pendulum at point 1
    x0=0;
    y0=0;
   %pendulum at point 2 
x1=-l*cos(results(i,1));
y1=-l*sin(results(i,1));
%plotting the cordinates
plot([y0 y1],[x0 x1],'marker','o')
%drawing a support for the pendulum using the line command
hold on
line([-2 2],[0 0])
hold off
%limiting the axis
axis([-2 2 -2 2])
%each frame developed by cordinates stored in m(ct)
pause(0.3)
M(ct)=getframe(gcf);
%ct is incremented at each itteration.
ct=ct+1;
end
%creating a movie
movie(M)
videofile=VideoWriter('pendulum_animation.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

figure(2)
%solve ode
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span,theta_0);
plot(t,results(:,2))
hold on
plot(t,results(:,1))
ylabel('plot')
xlabel('time')

function of the main code-

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

explaination-

1-all the variable are taken into considertion like

damping coefficient,d

gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball ,to calculate the

2-intial condition of the pendulum which has '0' displacement and with a angular velocity of 3 rad/sec, and the the time span

3-time span of 20 sec for a hundred value of time taken as t_span.

4.figure(1)

here we are animating the pendulum by apening for loop for length of t(span)

5-defining the values of the pendulum at state point 1 which is (0,0)

at state point 2 we calculate the (x1,y1) by taking the values of results w.r.t t_span. formula used

x1=-l*cos(results(i,1));
y1=-l*sin(results(i,1));

6-we plot the cordinate of two state point to get the pendulum position. we set the axis limit,we draw the support of the pendulum

by using the line command.

7-the frame of each itteration of theloop is stored in M(ct).

8-A movie is created using the command movie(M).

9-using the videoWriter command we name the file pendulum_animation,avi

10-we write the video file by combining the each frame into a single video file using the command writeVideo.

11-figure (2) shows the plot of the velocity and displacement w.r.t time.

the animation-

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

screenshot of error faced-

the error is due to missed semi colon in the ode function as it was showing error in the line 15 wich was a ode function.the error was fixed after cross cheking the each line in the ode_func.