solving second order odes
Introduction:
A simple pendulum consist of a metal mass held by a string to a stationary wall.The mass oscillates because of gravity and its inertia when lifted and then left to fall.This oscillation can be modeled using second order differential equation.
`b=0.05( Ns)/m`
`g=9.81 m/s^2`
`l=1m.`
`m=1kg.`
`theta=0 rad.`
`(thetadott)=3 (rad)/s`
Input equation
`(d theta1)/dt=theta2`
`(d theta 2)/dt=-(b/m)*theta2-(g/l)*sin(theta1)`
Main program
clear all
close all
clc
b=0.05
g=9.81
l=1
m=1
theta=[0;3]
t_span=linspace(0,20,150)
[t,results]=ode45(@(t,theta) odefunc(t,theta,b,g,l,m),t_span,theta)
figure(1)
hold on
plot(t,results(:,1))
plot(t,results(:,2))
xlabel('Time(s.)')
ylabel('Angular displacement and velocity(rad.)')
grid on
figure(2)
hold on
x0=0
y0=0
ct=1
for i=1:length(t)
x1=l*sin(results(i,1))
y1=-l*cos(results(i,1))
plot([x0,x1],[y0,y1],'color','r','linewidth',5)
hold on
plot([-1.4,1.4],[0,0],'linewidth',10)
hold on
plot(x1,y1,'o','markersize',20,'markerfacecolor','y')
axis([-1.4 1.4 -1.4 1.4])
grid on
hold off
pause(0.07)
M(ct)=getframe(gcf);
ct = ct+1;
end
movie(M)
videofile=VideoWriter('vib.avi','uncompressed avi')
open(videofile)
writeVideo(videofile,M)
close(videofile)
Main program explaination:
- Inputs are initialized
- Time array is generated using linspace command which starts at zero and ends at 20 with 150 points in between
- Then ode45 command is writted which outputs results which contain theta and theta dot with time and on the right side the main variable of equation is written first then the function which was made on another code with all its inputs and then the time span which it will integrate upon and then the initial conditions.
- then figure (1) plot is generated by plotting time with all the rows and first column for theta and and with all the rows and second column for theta dot
- Then figure 2 is the animation
- First x0 and y0 is defined
- Then counter is initiated
- Then for loop is constructed by using the length of the time array
- Then x1 and y1 is defined using geometric relation and uses the output angular displacement in column 1 and uses the counter ‘i’ for row counter
- Then the pendulum string is plotted using array brackets
- Then the wall is plotted using array brackets
- The the ball of the pendulum is plotted using a marker 0
- Then axis is defined
- then hold off was used to remove the plot in order to have only one plot on the screen after hold on was used to draw all the elements of the plot
- Then it is slowed down to 0.07
- Then the frame is captured
- Then counter is increased by 1
- Then movie is made in an uncompressed format.
Function program:
function [dtheta_dt]=odefunc(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
Function program explaination:
- Function output is defined which is theta dot array and then input variables are stated
- Then theta1 is defined as first element is an array
- Theta2 is defined as second element in an array
- Then equations ae declared
- Then the output is put in column matrix
programming error:
the main problem that happened was related to using hold on in order to plot all the elements together which when not used resulted in partial plot only and then using hold off to remove the whole plot in order to prevent multiple plots on the sceen . when not used this resulted in all plots plotted on top of each others.All the errors were then corrected.
Output plot:
Output animation: