Week 3 - Solving second order ODEs
A Mat lab program to simulate the transient behavior of a simple pendulum
In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a pendulum
The way the pendulum moves depends on Newton’s second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.
Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of these systems.
Your objective is 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.
Write a program in Mat lab\Octave that will simulate the pendulum motion, just like the one shown in the start of this challenge.
Use in program,
L=1 meter,
m=1 kg,
b=0.05.
g=9.81 m/s2.
Simulate the motion between 0-20 sec, for angular displacement=0,angular velocity=3 rad/sec at time t=0.
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
Main_Code:
clear all
close all
clc
%inputs
%Damping coefficient
b=0.05;
%gravity
g=9.81;
%length in meters
l=1;
% Mass in kg
m=1;
%intial condition
%one is angular velociity
theta_0=[0;5];
%time between o to 10 seconds of 250 values in between.
t_span=linspace(0,10,250)
[t,results]= ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);
%plotting
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel('Time (Sec)')
ylabel('Angular Displacement and Angular Velocity')
grid on
legend('Angular Displacement','Angular Velocity')
hold off
%plotting pendulum
figure
ct=1;
for i=1:length(t_span)
x_fixed=0;
y_fixed=0;
x_freeend=l*sin(results(i,1));
y_freeend=-l*cos(results(i,1));
plot([x_fixed x_freeend],[y_fixed y_freeend],[-0.5 0.5],[0 0],'linewidth',3)
hold on
plot(x_freeend,y_freeend,'marker','o','markerSize',45,'markerFacecolor','r')
axis([-1.5 1.5 -1.5 0.5])
grid on
pause(0.02)
hold off
M(ct)=getframe(gcf);
ct=ct+1;
end
movie(M)
videofile = VideoWriter('ode_pendulum.avi','uncompressed avi');
open(videofile)
writeVideo(videofile,M)
close(videofile)
Result:
