Week 3 - Solving second order ODEs

Program for second order differential equation (simple pendulum )

We take first input required to solve second differential equation for simple pendulum

b=0.05 (damping coefficient), g=9.81 (acceleration due to gravity), l=1 (length of string), m=1 (mass of pendulum)

Initial condition Is given as when t=0, angular displacement (theta) = 0 rad and angular velocity = 3 rad/sec for time span of 20 s.

theta_0 = [0;3]

We use linspace command to get 500 values in between 0 to 20 t_span = linspace(0,20,500)

We use ode45 function to solve differential equation, to convert higher order differential equation to first order differential equation we create new function as ode_func

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

to get the results

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

In above equation t and theta value given as the t_span and theta_0 as in initial condition

After getting results to plot graph of angular displacement and angular velocity vs time

figure(1) (for new plot window)

plot(t,results(:,1)) (: is for up to last value in row ),(1for first column)

hold on (for hold currant plot in same window )

plot(t,results(:,2)) (for second column results)

for simulation of pendulum we use angular displacement results to plot pendulum position

defining theta3 to store the results one

Use the for loop for getting different position of pendulum

for i=1:length(theta3)

THETA3 = theta3(i);

x0=0; (coordinates of initial point or fix point of pendulum)

y0=0;

x1=l*sin(THETA3); (coordinates of other end or free end of pendulum)

y1=-l*cos(THETA3);

figure(2) (for second window of plot)

plot([-0.5 0.5],[0 0],'color','g','linewidth',2); (plotting base line for pendulum)

hold on (hold in same window)

plot([x0 x1],[y0 y1],'color','r','linewidth',3); (plotting pendulum string)

hold on

plot(x1,y1,'o','markers',10,'markerfacecolor','b') (markers command for making shape of pendulum ball)

hold off

axis([-1.5 1.5 -1.5 0.5]) (shifting axis coordinates)

pause(0.001) (to give pause in between frames)

grid on

M(ct)= getframe(gcf); (to get frames for movie)

ct=ct+1; (count increment of frame)

link of video

https://youtu.be/ggrw8Fk66dI

close all
clear all
clc

% Program for second order differential equation
%Inputs
b=0.05;
g=9.81;
l=1;
m=1;

% initial condition 
theta_0=[0;3];

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

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

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('Angular velocity & Angular displacement')

theta3= results(:,1);
ct= 1;
for i=1:length(theta3)
    THETA3 = theta3(i);
    
    x0=0;
    y0=0;
    
    x1=l*sin(THETA3);
    y1=-l*cos(THETA3);
    
    figure(2)
    plot([-0.5 0.5],[0 0],'color','g','linewidth',2);
    hold on
    plot([x0 x1],[y0 y1],'color','r','linewidth',3);
    hold on
    plot(x1,y1,'o','markers',10,'markerfacecolor','b')
    hold off
    axis([-1.5 1.5 -1.5 0.5])
    pause(0.001)
    grid on
    M(ct)= getframe(gcf);
    ct=ct+1;
end
movie(M)
Videofile = VideoWriter('Pendulum_ODE.avi','Uncompressed AVI');
open(Videofile)
writeVideo(Videofile,M)
close (Videofile)

ODE 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