Week 3 - Solving second order ODEs
objective = To write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping
In the above equation,
g = gravity in m/s2,
L = length of the pendulum in m,
m = mass of the ball in kg,
b=damping coefficient.
Input Values
L=1 metre,
m=0.1 kg,
b=0.05.
g=9.81 m/s2.
Convert 2nd degree ODE order into first degree ODE
Ɵ = Ɵ1
dƟ/dt = dƟ1/dt = Ɵ2
dƟ/dt+b*Ɵ/m+g*sinƟ1/L
Create an 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
Program of a simple pendulum -
clear all
close all
clc
b = 0.05; % damping co-efficient
m = 0.1; % mass - kg
g = 9.81; % gravity - m/(s^2)
l = 1; % length - meter (m)
% 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)
i = 1
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('time')
ylabel('plot')
for i = l:length(t)
x = (sin(results(i,1)))
y = (-cos(results(i,1)))
figure(2)
plot([0, x],[0, y],'-o')
axis ([-1 1 -2 0])
text(0, -1.6, ['time', num2str(t(i,1))])
M(i) = getframe(gcf);
pause(0.1)
end
movie(M)
videofile=VideoWriter("simple_pendulam.avi","uncompressed AVI");
open(videofile)
writeVideo(videofile,M)
close(videofile)
Result -
plot-
Animation video -
https://drive.google.com/file/d/1Qr2jZh0RFKC31GUqvrIqXtIRkV-vtiVI/view?usp=sharing