Week 3 - Solving second order ODEs

https://drive.google.com/file/d/14oS4FxA9mFT-_0E0jq_fHbWchgcFoq6Y/view?usp=sharing

The above link given for animation of simple pendulum

Explanation for each step of program given at the bottom

clear all
close all
clc

b=0.05;
g=9.81;
l=1;
m=1;

% initial condition
theta_0 = [0;3];

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

% solve ode

[t, results] = ode45(@(t, theta) func(t, theta, b,g,l,m), t_span, theta_0) 
y =  results(:,1);
plot(t,results(:,1));
hold on;
plot(t,results(:,2));
xlabel('time');
ylabel('plot')
ct=1;
 
for i=1:length(y)
    x0=0;
    y0=0;
    x1 = l*sin(y(i));
    y1= -l*cos(y(i));
    plot([-1.5 1.5],[0 0])
    hold on
    plot([x0 x1],[y0 y1],'linewidth',3,'color','r')
    
    plot(x1,y1,'marker','o','markersize',20,'markerfacecolor','g')
    
    title({'simple harmonic motion';'of a pendulum'},'color','r')
    axis([-1.5 1.5 -1.5 1.5])
    pause(0.01);
    M(ct)=getframe(gcf);
    hold off
    ct=ct+1;   
end

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

close all
clc

b=0.05;
g=9.81;
l=1;
m=1;

% initial condition (we are initializing the condition here)

theta_0=[0;3]; (Here zero displacement and 3rad/sec is given as initial angular velocity)

% time points
t_span=linspace(0,20,20); (Here we are giving the time limit between 0sec and 20 sec for the integration and 20 intermediate steps also giving)

% solve ode

[t,result]=ode45(@(t,theta)func(t,theta,b,g,l,m),t_span,theta_0);

(Here we are solving the ode by using ode45 command which is inbuilt in MATLAB and for that 'func 'is using and we are giving the inputs accordingly where 't_span' is the time array and 'theta_0' is the initial condition. For 'func' call we are writing seperate program)

y=result(:,1); (Here we are taking the 1st coloumn of results (ie angular displacement) and assigning that to 'y')

plot(t,result(:,1)); (Here we are plotting time versus 1st column of result array (anguar displacement) vs time)
hold on; (It is giving otherwise only latest plot only will show)
plot(t,result(:,2)); (Here we are plotting time versus 2nd column of result array (anguar velocity) vs time)
xlabel('time'); (Here we are giving name for the xaxis of the plot as 'time')
ylabel('plot') (Here we are giving name for the yaxis of the plot as 'plot')

ct=1; (Here we are initializing the frame count)
for i=1:length(result(:,1)); (Here we are opening the for loop for finding different theta positions by giving 1 to 20 values for 'i')
x0=0; (Giving starting point of pendulum)
y0=0; (Giving starting point of pendulum)

x1=l*sin(y(i)); (Finding different values for x1)
y1=-l*cos(y(i)); (Finding different values for y1)
plot([x0 x1],[y0 y1],'linewidth',3,'marker','o','color','r') (Plotting pendulum leg and 'o ring' a the end of pendulum)
hold on
plot(x1,y1,'marker','o','markersize',20,'markerfacecolor','g') (Here we are giving the size of that end point)
hold on
title({'simple harmonic motion';'of a pendulum'},'color','r') (Here we are giving tiltle)
axis([-1.5 1.5 -1.5 1.5]) (Here we are fixing the axes)
pause(0.01); (Giving the pausing time for each frame)
M(ct)=getframe(gcf); (Storing all frames)
ct=ct+1; (Increasing the frames)

end (Closing the for loop)

% creating animation for S.H.M (Here we are creating the animation video)
movie(M)
videofile=videoWriter('pendulum.avi','Uncompressed AVI'); (Wring video by using videowriter and save that into videofile)
open(videofile); (Opening that videofile)
writeVideo(videofile,M);
close(videofile) (Closing that videofile)

Separe program written for 'func' which is given below:

function[dtheta_dt]= func(t,theta,b,g,l,m)

(Here 'dtheta_dt' is the function handle , 'ode_func' is the function name and t,theta,b,g,l and m are the corresponding inputs)

theta1=theta(1); (Here we are defining the 'theta1')
theta2=theta(2); (Here we are defining the 'theta2')

dtheta1_dt=theta2; (Here we are giving dtheta1_dt=theta2)
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1); (Here we are defining dtheta2_dt)
dtheta_dt=[dtheta1_dt;dtheta2_dt]; (Here we are combining the dtheta1_dt and dtheta2_dt)
end (Ending the function)

while doing this program I got some errors in command window , that errors I am attaching with this.

For solving error 1, I given function handle name as [dtheta_dt]. Before getting this error I didn't given that function handle name.

For solving error 2, I given (t,theta) while writing that ode45 command line. Before getting this error I wrote (t,result) instead of writing (t,theta)

After run this program I got a plot which shows the variation of angular displacement and angular velocity with time and that plot also I am attaching along with this.