Week 3 - Solving second order ODEs

%simple pendulum
clear all
close all
clc
 
%inputs
b = 0.05;
g = 9.81;
l = 1;
m = 1;
 
% initial condition
theta_0 = [0;3];
 
% time points
t_span = linspace(0,20,500);
 
% solve ODE
[t, results] = ode45(@(t,theta) ode_func(t,theta,b,g,l,m),t_span, theta_0);
 
% ploting
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('time in sec')
xlabel('velocity/displacement')
 
% Iteration using for loop to create animation of simple pendulum
ct=1;
theta_1 = results(:,1);
 
for i= 1:length(results(:,1))
    
    % initial co-ordinates
    x0= 0;
    y0= 0;
    theta_2 = theta_1(i);
    
    % final co-ordinates
    x1 = 1*sin(theta_2);
    y1 = -1*cos(theta_2);
    
    % ploting for animation
    figure(2)
    % base of pendulum
    plot([x0 x1],[y0 y1],'linewidth',3,'color','black');
    hold on
    plot(x1,y1,'o','markers',15 , 'MarkerFaceColor','green');
    hold on
    plot([-3 3],[0 0],'linewidth',3,'color','red'); 
    hold off
    axis([-3 3 -2 2])
    grid off;
    pause(0.03);
    M(ct) = getframe(gcf);
    ct = ct+1;
end
% creating movie from captured frames
 movie(M)
 videofile = VideoWriter('simple_pendulum.avi','uncompressed AVI');
 open(videofile);
 writeVideo(videofile,M);
 close(videofile);

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_dt; dtheta2_dt];
end

OUTPUT:-

link for simulation video of simple pendulum:-

https://youtu.be/2HKdjA_ByJE

explanation of code:-

1. Now we are solving the ODE as mentioned above and we split that ode into two single order ODE’s one for velocity and one for displacement.
2. So, when we need to solve the first ODE, we need to provide the initial conditions for the displacement and velocities.
3. Since in MATLAB “ode45” is an inbuilt solver, which is going to solve the ODE system.
4. Our ODE depends upon the time, theta and other parameters i.e.,

• Damping coefficient (b)
• Acceleration due to gravity (g)
• Length (l)
• Mass of bob (m)

5. For making bob I use the marker command. And keeping the one end of the pendulum fixed so that we can oscillate our pendulum
6. Initially used commands - 'clear all' used to clear all the assigned values to the variables.

   'close all' - force deletes all figures

   'clc' - clears the command window

7. Initial conditions which are displacement and velocity are given in the variable 'theta_0'. These conditions are further used to solve the ODE

8. The next step includes giving the time points for the ODE. These points are used for integrating while solving the ODE. The 'linspace' command is used to divide the time interval into required parts.

9. The function ' ode_fun ' is callled for solving the ODE using the 'ode45' command. Using this function we get Time, Displacement and Velocity values. A second order ODE is converted in first order ODE. The general syntax of the 'ode45' is as follows for this code is as follows:-

   [t,results] = ode45(@(t,theta) ode_func (t,theta,b,g,l,m),t_span, theta_0)

10. Plotting the 'Velocity & displacement vs Time' graph. The 'plot' and the 'hold' commands are used to plot the graph.

11. To create an animation for the pendulum a conter variable 'ct' which is equaled to 1 is created. Then a for loop has to be created from 1 to the length of time span and results array generated earlier.To do this the 'length' command is used. The following main code was used while plotting

12. The 'getframe' command is used to create a movie. For this the counter variable 'ct' is used to generate the movie array.

13. The created movie is used to create a video file using the 'VideoWriter' command. The 'open' and 'close' commands are used to create the open and close the file.Then 'writeVideo' is used to write the video into the file.The file is saved in the '.avi' format.

    ERRORS:-

    I wrote 'o' instead of '0' in ode45 line so I got an error. I corrected it by replacing 'o' with '0 '