Week 3 - Solving second order ODEs

function[dtheta_dt]=ode_func(t,theta,b,g,l,mass)

theta1=theta(1);
theta2=theta(2);
dtheta1_dt=theta2;
dtheta2_dt=-(b/mass)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

% Second order ODE function
clear all
close all
clc

%inputs
b=0.05;%damping coefficient
g=9.81;%acceleration due to gravity
l=1;%length of the string
mass=0.1;%mass of the pendulum

% Initial condition
theta_0=[0;5];

% time points
t_span=linspace(0,10,800);

% Solve ODE
[t,y]=ode45(@(t,theta)ode_func(t,theta,b,g,l,mass),t_span,theta_0);

figure(1)
plot(t,y(:,1))
hold on
plot(t,y(:,2))

ylabel('angular displacement and angular velocity')
xlabel('time')
legend('angular displacement','angular velocity')

%% animation
ct=1;
for i=1:length(t_span)
    x0=0;
    y0=0;

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

    %pendulum motion simulation

    figure(2)
    plot([x0 x1],[y0 y1],'linewidth',2)
    hold on
    axis([-2 2 -1.5 0.5])
    viscircles([0 0],0.05);% to create a rigid hook
    viscircles([x1 y1],0.20);% to create a pendulum ball
    
    title('motion of a simple pendulum')
   

    Mov(ct)= getframe(gcf);
    ct=ct+1;
    hold off
end

%animation of motion of simple pendulum
movie(Mov)
videofile=VideoWriter('simple_pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,Mov)
close(videofile)

YOUTUBE LINK FOR THE VIDEO:

https://youtu.be/nShlacK142I

EXPLANATION OF CODE:

1.start the code with clear all,close all,clc

2.the inputs are given as follow

%inputs

b=0.05;%damping coefficient

g=9.81;%acceleration due to gravity

l=1;%length of the string

mass=0.1;%mass of the pendulum

with x0 and y0 being 0,0 and x1,y1 being

3.x1=1*sin(y(i,1));

y1=-1*cos(y(i,1));

4.the initial condition of theta starts from 0 to 5.

5.the timespan is given from 0 to 10 with 800 divisions using linspace.

6.before solving the 2nd order ode we have to declare the function in ode45

7.the function is used to formulate the 2nd order ode formula such as

function[dtheta_dt]=ode_func(t,theta,b,g,l,mass)

theta1=theta(1);

theta2=theta(2);

dtheta1_dt=theta2;

dtheta2_dt=-(b/mass)*theta2-(g/l)*sin(theta1);

dtheta_dt=[dtheta1_dt;dtheta2_dt];

end

8.we know that

d2theta/dt2+b/m(dtheta/dt)+g/l(sin(theta))=0

b=length

m=mass

g=gravity

l=length

theta=theta(1)

dtheta/dt=dtheta1/dt=theta2

d2theta/dt2=d2theta1/dt2=d/dt[dtheta1/dt]=d/dt(theta2)

9.d2theta/dt2+b/m(dtheta/dt)+g/l(sintheta)=0

d2theta1/dt2+b/m(dtheta1/dt)+g/l(sintheta1)=0

dtheta2/dt+b/m(theta2)+g/l(sintheta1)=0

dtheta1/dt=theta2

dtheta2/dt=(-b/m)-(g/l)sintheta1

10. the above formula was formulated in function and was declared in the ode45.

the ode45 function is used to solve the two input arguments such as t and y.

by declaring the function as follows

[t,y]=ode45(@(t,theta)ode_func(t,theta,b,g,l,mass),t_span,theta_0);

11.we now go to plot the ode and animate it

by using the below command

figure(1)

plot(t,y(:,1))

hold on

plot(t,y(:,2))

ylabel('angular displacement and angular velocity')

xlabel('time')

legend('angular displacement','angular velocity')

12.%% animation

ct=1;

for i=1:length(t_span)

x0=0;

y0=0;

x1=1*sin(y(i,1));

y1=-1*cos(y(i,1));

13.%pendulum motion simulation

figure(2)

plot([x0 x1],[y0 y1],'linewidth',2)

hold on

axis([-2 2 -1.5 0.5])

the axis starts from -2 to 0.5

viscircles([0 0],0.05);% to create a rigid hook

viscircles([x1 y1],0.20);% to create a pendulum ball

14.viscircles is used for giving definite radii to a circle and shape

title('motion of a simple pendulum')

15.Mov(ct)= getframe(gcf);

ct=ct+1;

hold off

end

%animation of motion of simple pendulum

movie(Mov)

videofile=VideoWriter('simple_pendulum.avi','Uncompressed AVI');

open(videofile)

writeVideo(videofile,Mov)

close(videofile)

errors i faced:

1.the string was not moving towards x axis but it was only moving in y axis

resolving the error:

1.the points i declared was not correct

such as

plot([x0 y1],[y0 y1],'linewidth',2)

i declared y1 instead of x1

corrected one-

plot([x0 x1],[y0 y1],'linewidth',2)