Week 3 - Solving second order ODEs

1. AIM: To write a program to solve a second order ordinary differential equation, simulate and animate a pendulum using MATLAB.

2. GOVERNING EQUATIONS:

Second order ordinary differential equation used to simulate the motion of a pendulum with large angular displacement or amplitude:

`(d^2 θ)/(dt^2 )+b/m(dθ)/dt+g/L sin⁡θ=0`

b= damping coefficient

m=mass

g=gravity

l=length

Conversion to first order ODE:

Therefore:

3. OBJECTIVES OF THE PROJECT:

1. To write a function for second order ODE.
2. Write a program to simulate the movement of the pendulum and call the function.
3. Animate the movement of the pendulum

4. BODY OF THE CONTENT:

clear
clc
close all 
% declaring constants
b = 0.05;
g = 9.81;
l = 1;
m = 1;
 
% initial condition for displacement and velocity
y0 = [0;3];
 
% time points for the duration of 20s
t_span = linspace(0, 20, 100);
 
%solve ODE
[t, y] = ode45(@(t, theta) ode_func(t, theta, b, g, l, m),...
              t_span,y0);
 
% plot angular velocity and displacement with time          
figure (1)
plot(t,y(:,1),'r')
hold on
plot(t, y(:,2),'g')
xlabel('Time (s)')
ylabel('\theta/\omega')
 
z = y(:,1);  %z stores all values from column 1
 
ct = 1; % start counter at 1
 
for i = 1:length(z)  % repeat for loop upto the length of z which is 100
    I = z(i);
    x0 = 0;
    y0 = 0;
    x1 = l*sin(I);  % x coordinate of free end of pendulum
    y1 = -l*cos(I); % y coordinate of free end of pendulum    
    figure (2)
    % plot fixed horizontal line
    plot([-0.5,0.5],[0,0], 'linewidth', 5, 'color', [0.9290 0.6940 0.1250]) 
    line([x0 x1],[y0 y1])
    hold on;
    % plot free end of pendulum
    plot(x1,y1,'o','MarkerSize',10,'MarkerEdgeColor','red',...
    'MarkerFaceColor',[1 .6 .6]);
    grid on ;
    axis([-1.5 1.5 -1.5 0.5]);
    pause(0.1);
    hold off;
    M(ct) = getframe(gcf);
    ct = ct + 1;
end
    
% creating movie similar to robot arm challenge
movie(M)
videofile = VideoWriter('oscillating_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Function:

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

1. First the constants are decleared.
2. Initial values of displacement and velocity are given in a column matrix form.
3. Time duration from 0 to 20 sec. is considered with 100 equal steps.
4. Then the function is solved by writing it in the form [t,y] = ode45(odefun,tspan,y₀), where t and y are the output arguments. tspam is the time limin and y₀ are the initial conditions. odefun is the function that is to be solved.
5. ode45 only works with functions that use two input arguments. However, we can pass in extra parameters by defining them outside the function and passing them in while specifying the function handle.
6. Here parameters other than theta and t are defined in a separate function called ode_func.
7. The ode45 solves the function to return velocity and displacement as output for 100 values of time between 0 to 20 s.
8. This velocity and displacement is plotted against time in figure 1.
9. Next we need to plot 100 positions of the pendulum for corresponding values of y.
10. For loop is used for doing so.
11. The animation is made using exact same procedure used in Robot Arm challenge.
12. For the function ode_func, the second order ODE is first converted into first order ODE and ode_func represents this system of equations as a function that accepts six input arguments: t, theta, b, g, l, and m.

6. RESULTS:

https://youtu.be/SyDx-0yXdEE

6. CONCLUSION:

Ode45 is very handy function to calculate linear ordinary differential equations.

7. REFRENCES:

https://www.acs.psu.edu/drussell/Demos/Pendulum/Pendulum.html

https://in.mathworks.com/help/matlab/ref/ode45.html

https://in.mathworks.com/help/matlab/math/parameterizing-functions.html

https://in.mathworks.com/help/matlab/math/choose-an-ode-solver.html

https://in.mathworks.com/help/matlab/math/summary-of-ode-options.html