Week 4.1 - Solving second order ODEs

Aim- WAP to solve the given ODE function for simple pendulum with damping.

Governing Equation-

dtheta2/dt = -(b/m)*theta2 - (g/l)*sin(theta1)

Given -

g = 9.81 (gravity)

L = 1 (length of the pendulum)

m = 1 (mass of the ball)

b = 0.05 (damping coeffient)

Body of Content-

1. Created a funtion file

2. used for loop

3.used plot command

Code-

%%second order ode %%
%% main code%%
clear all
close all
clc

%% inputs%%

%%length of the pendulum%%
l=1;
%%mass%%
m=1;
%%damping coefficient%%
b=0.05;
%%gravity%%
g=9.81;

%initial condition%%
theta_0 = [0;3];

%%time point%%
t_span = linspace(0,20,100);

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

figure(1)
title('parameter variation w.r.t. time plot')
plot(t,results(:,1))
hold on 
plot(t,results(:,2))
ylabel('plot')
xlabel('time')

legend('displacement','velocity')
grid on
ct=1;
angle= results(:,1);

for i=1: length(t_span)
    Angle =angle(i);
    figure(2)

    x0=0;
    y0=0;
    y1=-1*cos(Angle);
    x1=1*sin(Angle);

    %%plots%%

    plot([x0 x1],[y0 y1],'linewidth',4)
    hold on
    plot(x1,y1,'.','Markers',50)
    hold on 
    plot([-5,5],[0,0],'LineWidth',5)
    hold off
    axis([-9,9,-4,4])
    pause(0.01)
    M(ct) = getframe(gcf);

    ct = ct+1;
end

movie(M)
videofile = VideoWriter('pendulum anime.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

%%ode function%%
%% code for function%%

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

theta1 = theta(1);
theta2 = theta(2);
dtheta1_dt = theta(2);
dtheta2_dt = -(b/m)*theta2 - (g/l)*sin(theta1);
dtheta_dt = [dtheta1_dt;dtheta2_dt];
end

Conclusion : Written a program and created animation for simple pendulum.

Link for the animation-

https://drive.google.com/file/d/1QnsE9CS68WOw5bqfND_8RgOYwqBu8v_D/view?usp=sharing