Week 3 - Solving second order ODEs

% Program to simulate the pendulum motion

clear all
close all
clc

%inputs

b = 0.05;               % Damping coefficient
g = 9.81;               % gravity constant
l = 1;                % Length of pendulum
m = 1;                % Mass of ball

theta_0 = [0;3];        % two initial condition,ie, displacement and velocity for two ODEs which depend on them individually

t_span = linspace(0,20,400);        %time array which goes from 0 to 20 (is divided into 300 parts)

%solve ODE
[t, results] = ode45(@(t, theta) ode_func(t, theta, b, g, l, m), t_span, theta_0);      % syntax for defining inputs and ode function,ode_func, with inbuilt solver,'ode45', in Matlab

for (i=1:length(t))

    P = l * [sin(results(i,1))  -cos(results(i,1))];            % coordinates of pendulum ball taking origin as the hinge point of string
    
    pend_ball = plot(P(1), P(2),'marker','.' ,'MarkerSize', 50, "Color",'b');    % plotting the pendulum ball 

    pend_string = plot([0 P(1)], [0 P(2)], 'linewidth', 3, 'color', 'k');        % plotting the string of the pendulum

    hold on                     % for animation in multiple 'plot' since plot command displays fig of most recent plot execution in multiple plot commands 
    pause(0.02)                 % pausing execution to get all the frames' execution in same window
    axis ([-1 1 -1.1 0.2]);     % changing the origin to properly accomodate the pendulum
    grid on                     % for making grid in the plot fig

    M(i) = getframe(gcf);        % M array getting each frame for video

    if (i<length(t))            % for deleting the plots after plotting it in figure before the next cycle
        delete(pend_string);
        delete(pend_ball);
    end
end

movie(M)                             %codes for saving the animation in a video file
videofile = VideoWriter('pendulum motion.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

Below is the code of the function 'ode_func' which is used for solving second order ode by changing it into first order ode.

function [dtheta_dt] = ode_func(t,theta,b,g,l,m)

    theta1 = theta(1)       
    theta2 = theta(2)
    dtheta1_dt = theta2;                                % Changing the 2nd ODE into two 1st ODEs for solving it
    dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);     
    dtheta_dt = [dtheta1_dt; dtheta2_dt];               % Output as two rows of which first is distance vector and second is velocity vector
end

In this program, the real motion of a pendulum with damping is created by solving its ordinary differential equation(ode). To solve a 2nd order ode in matlab, we have to change it into a 1st order ode which is done in the function 'ode_func'. Most of the explanation of code is done in the code itself.

To solve an ode in matlab we use its in-built ode solver ode45. ode45's first input defines the main inputs for the ode function ,i.e., time 't' and theta. The second input defines the function 'ode_func'. The third input gives the initial condition, integration of functions and time array for the length of the integration. The fourth input defines the initial condition at the start of the time array.