Week 3 - Solving second order ODEs

1. AIM:

To solve a second order ordinary differential equation in matlab.

2. GOVERNING EQUATIONS:

Variables:

- θ is the angle made by rod with vertical

- b is damping coefficient

- m and L are mass of the bob and length of the rod / string

b = 0.05(Ns)/m

l=1m

m=1kg

θ=0rad.

(θt')=3rad/s

Input equation:

dθ1/dt=θ2

dθ2dt=−(b/m)⋅θ2−(g/l)⋅sin(θ1)

4. Main program:

% Inputs 
 % Damping coefficient
 b = 0.05;
 % Acceleraton due to gravity - m/s^2
 g = 9.81;
 % Mass of the bob - kg
 m = 1;
 %Length of the string - m
 l = 1;

% Initial condition 
dis_angvelo = [0;3];

% Timepoints
t_span = linspace(0,20,500);

% Solving ode 
[t,results] = ode45(@(t,theta) ode_function(t,theta,b,g,l,m),t_span,dis_angvelo);

% Plotting 
figure(1)
plot(t,results(:,1),'linewidth',1)
hold on 
plot(t,results(:,2),'linewidth',1)
axis([0 25 -3 3])
xlabel('Time - sec')
ylabel('Angular Velocity & Displacement')
legend('change in displacement','change in angular velocity')
hold off

figure(2)
plot(results(:,1),results(:,2))
xlabel('Position - rad')
ylabel('Angular Velocity - rad/sec')


theta2 = results(:,1);
ct = 1;
for i = 1:length(theta2);
    theta3 = theta2(i);
    
    x0 = 0;
    y0 = 0;
    
    x1 = l*sin(theta3);
    y1 = -l*cos(theta3);
   
    figure(3)
    plot([-1.1 1.1],[0 0],'linewidth',2);
    hold on 
    plot([x0 x1],[y0 y1],'color','r','linewidth',2)
    hold on 
    plot(x1,y1,'o','markers',20,'markerfacecolor','b')
    hold off 
    grid on 
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.003)
    M(ct) = getframe(gcf);
    ct = ct+1;
end
movie(M)
videofile = VideoWriter('pendulum.avi','uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)

5. Main program explaination:

• Inputs are initialized
• Time array is generated using linspace command which starts at zero and ends at 20 with 500 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.003
• Then the frame is captured
• Then counter is increased by 1
• Then movie is made in an uncompressed format.

6. 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.

• The string or the rod of the pendulum was not visible. Fixed it by putting the marker dot in the initial plot itself, and not in another.

7. PLOTS FOR THE SIMPLE PENDULUM :

Plot for Time and AngularVelocity, Displacement.

8. Animation for simple pendulum :-

https://youtu.be/XOtZuc9WddY