Week 3 - Solving second order ODEs

Function for 2nd order ODE :-

% function for 2nd order ode 

function [dtheta_dt] = ode_simple_pendulum(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

Program to solve 2nd order ODE taking simple pendulum as example :-

%--------------------SIMPLE PENDULUM-------------------%
% Program to solve a 2nd order ODE (d^2theta/dt^2)+(b/m)*(dtheta/dt)+(g/L)*sin(theta)=0
    

clear all
close all
clc


b = 0.5;   % b is the damping factor
g = 9.81;  % g is the acceleration due to gravity
L = 1;     % L is the length of the rod
m = 1;     % m is the mass attached to the rod

% Initial conditions
S_V = [0;3];          % S is the angular dispalcement, V is the angular velocity

% timepoints
t_span = linspace(0,20,500);  % taking time period of 0 to 20 sec for 1000 points

% solve ode
[t, results] = ode45(@(t,theta) ode_simple_pendulum(t, theta, b, g, L, m), t_span,S_V);


% Plot for angular displacement and angular velocity w.r.t time

figure(1)
plot(t,results(:,1))
hold on
grid on
plot(t, results(:,2))
xlabel('time')
ylabel('plot')
title('projection of angular velocity and angular dispalcement w.r.t time for simple pendulum')
legend('displacement', 'velocity')


% Creating an animation to demonstrate simple pendulum using for loop and get frame

theta2 = results(:,1);
ct = 1;

for i = 1:length(theta2)
    THETA_1 = theta2(i);

    X0 = 0;
    Y0 = 0;

    X1 = L*sin(THETA_1);
    Y1 = -L*cos(THETA_1);
figure(2)
    plot([-0.2 0.2], [0 0], 'color', 'r','linewidth',3) % plot for fixed end
    hold on
    plot([X0 X1], [Y0 Y1], 'color', 'b', 'linewidth',1.5) % plot for rod connecting fixed end and mass
    hold on
    plot(X1,Y1,'O','markers',10,'markerfacecolor', 'y') % plot for a spherical mass 
    hold off
    grid on
    axis([-2 2 -1.2 0.1]) % taking axis dimensions from xaxis -2 to 2 and y axis from -1.2 to 0.1
    pause(0.003)
    M(ct) = getframe(gcf);
    ct = ct+1;
end

% Creating a videofile
movie(M)
videofile = VideoWriter('Simple_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Plot :-

Animation link :-

https://youtu.be/0DhXxlnW6mA

Explanation of Code :-

Step_1 :-

-> To solve a 2nd order ODE, we are creating a function.

-> In this function we are converting 2nd order ODE into 1st order ODE to get two results and assigning them into an array dtheta_dt.

-> IN order to covert we taking theta _1 = theta(1), theta_2 as theta(2).

Step_2 :-

-> Now we are writing the main code to solve the 2nd order ode.

-> In this step we are giving the inputs and initial conditions provided in the question.

-> Also taking the time span of pendulum oscillations up to 20sec for 500 values.

Step_3 :-

-> Writing ode45 syntax

[TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0)

-> ODEFUN is a function handle. For a scalar T and a vector Y, ODEFUN(T,Y) must return a column vector corresponding to f(t,y).

-> TSPAN is a two-element vector [T0 TFINAL] or a vector with several time points [T0 T1 ... TFINAL].

-> YO is a column vector of initial conditions, one for each equation.

Step_4 :-

-> Plotting the results in a single figure using hold on command. And labeling the plot.

Step_5 :-

-> Creating an animation to demonstrate simple pendulum using for loop and get frame.

-> taking theta_2 = results(:,1), giving ct = 1;

-> Creating a for loop and taking i = length(theta2).

-> Assigning theta2(i) to THETA1.

-> giving the initial values x0 = 0; y0 = 0;

-> and taking the end point x1 = L*sin(THETA1); y1 = -L*cos(THETA1);

Step_6:-

-> Now we are creating a plot to get the frames.

-> firstly, we are plotting the fixed points and giving hold on command.

-> for the second one we are plotting the length of the rod.

-> for the third one we are creating a plot to make a circular anime using marker command which is considered as weight connected to the rod.

-> Now setting the axis using axis command and giving a pause limit of 0.003.

-> Creating the movie frame using getframe command and giving ct = ct+1, ending the for loop.

Step_7 :-

-> After the plot is converted into frames now, we are creating it into a movie(M).

-> Creating a video file by writing it as simple_pendulum.avi which is uncompressed avi.

-> Hence video file is created and we can open it as given in the youtube link.

Error's found :-

Explanation of error :-

-> In the above fig we can see the animation having lot of frames which are fixed while oscillating.

-> this error is found due to misplacing of hold on as marked with red ink as shown in above fig.

-> This this way of presentation will be not good so we need to place the hold on command at write place to overcome those errors.