Week 3 - Solving second order ODEs
%ode function program
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;dtheta2_dt]
end
%program for simple pendulum
clear all
close all
clc
%input variables
b=0.05
g=9.81
l=1
m=1
theta_0=[0;3]
t_span=linspace(0,20,50)
%ode 45 function
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0)
%plotting of angular velocity and displacement
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
grid on
title('motion of pendulum')
xlabel('time in sec')
ylabel('angular velocity/displacement')
legend('displacement','angular velocity')
angle=results(:,1)
%count command
ct=1
for i=1:length(angle)
omega=angle(i)
figure(2)
title('motion of pendulum')
%fixed coordinates
x0=0
y0=0
%moving coordinates
x1=l*cos((3*pi/2)+omega)
y1=l*sin((3*pi/2)+omega)
%plotting
plot([x0 x1],[y0 y1],'linewidth',3,'color','g')
hold on
plot(x1,y1,'o','markers',15,'MarkerFacecolor','k')
hold on
plot([-3 3],[0 0],'linewidth',2,'color','b')
hold off
axis([-3 3 -2 2])
grid on
pause(0.1)
M(ct)=getframe(gcf)
ct=ct+1
end
%animation of simple pendulum
movie(M)
videofile=VideoWriter('simple_pendulum.avi','uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)
Explanation:
step 1: First create a function program and save it with the file name "ode_func'.
step 2: start the program for simple pendulum.
step 3: give the input variables.
step 3: Give the ode 45 function.
step 4: Give the inputs for plotting angular velocity and displacement.
step 5: Give the count commant for iteration.
step 6: Give the fixed coordinates
step 7: Give the moving coordinates
step 8: Now, give inputs for plotting.
step 9: Don't forget to put 'end' after this.
step 10: Give the inputs for animation of simple pendulum.