Week 3 - Solving second order ODEs

Function Code

%function to be called for pendulum animation
function [dtheta_dt] = ode_func(t,theta,c,g,l,m)
theta1=theta(1)
theta2=theta(2)
dtheta1_dt=theta2;
dtheta2_dt= -(c/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt= [dtheta1_dt;dtheta2_dt];
end

Main Code

close all
clear all
clc

c=0.05;
g=9.81;
l=1;
m=0.1;

%initial condition
theta_0 = [0;5];

%time points
t_span = linspace(0,20,500)

%solveODE
[t,results] = ode45(@(t, theta) ode_func(t,theta,c,g,l,m),t_span,theta_0);
plot(t, results(:,1))
hold on
plot(t,results(:,2))
% subplot(5,4,[15,16,19,20]);
% plot(results(:,1), results(:,2));
% xlabel('radians');
% ylabel('angular velocity');

 ct=1;
 for i= 1:length(t_span)
    THETA1(i)=results(i,1);
    x0=0;
    y0=0;
    x1= l*sin(THETA1(i));
    y1=-l*cos(THETA1(i));
    % subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);
    plot([x0 x1],[y0 y1],'color','k');
    hold on
    line([-1 1], [0 0],'color','g');
    plot(x1,y1,'o', 'markers',10,'markerfacecolor','r');
    xlim([-1.1, 1.1])
    ylim([-1.1,0.5])
    xlabel(' x-coordinates (m)')
    ylabel(' y-coordinates (m)')
    title("Pendulum Animation")
    grid on;
    axis([-1.5,1.5,-1.5,0.5])
    pause(0.05);
    hold off
    q(ct)= getframe(gcf);
    ct=ct+1;
end

movie(q)
videofile = VideoWriter( 'Pendulum_animation.avi', 'Uncompressed AVI');
open(videofile)
writeVideo(videofile,q)
close(videofile)

Coding logic - a function is used to break down the second-order ODE into two single order ODEs. In another program, the function is called inside 'ode45' for solving the ODEs. Then a for loop is used to iterate the horizontal and vertical displacement of the pendulum into angular displacement and put it in motion by plotting it.

For Loop Explanation - the loop is iterated till the length of the time points. A starting and ending x and y value is calculated and iterated till 'theta' values for the angular displacement. Now the same is plotted to get the animated effect of the pendulum. The marker command is used to create the bob of the pendulum.

Syntax ode45: [t,x] = ode45(odefun, tspan, x0)

Pendulum Animation link - https://drive.google.com/file/d/1CMmuGVVRQQ8mJIbjunAOOpr81IytljTx/view?usp=sharing