Week 3 - Solving second order ODEs
Here We can simulate the transient behaviour of a simple pendulum and to create an animation of it's motion.
Below is the figure of simple pendulum :
Below is the Second order differential equation for simple pendulum:
Input Values:
b=0.05 (Damping coefficient)
g=9.81 (gravity in m/s^2)
l=1 (Length of pendulum in meter)
m=1 (Mass in Kg)
In intial condition Lets take,
Assume,
theta_0=[0;4];
Angular Dislplacement is 0 and Angular Velocity is 4 and
t_span=linspace(0,20,300);
Time span is 0 to 20 sec, with 300 intervals.
Now solve ODE,
In matlab ode45 is the inubuilt solver for ODE,
Here ODE is the function of t and theta , t and theta are dependent and independent variables.
Below is the syntax for ODE,
[t,result]=ode45(@(t,theta) ode_function(b,g,l,m,t,theta),t_span,theta_0)
We have called function of name ode_function into it and give a name to ode_function as [dtheta_dt] First Simplify,
Assume,
theta=theta1................(1)
So, dtheta_dt=dtheta1_dt
dtheta1_dt=theta2.................. (2)
Then, d^2theta_dt^2=d^2theta1_dt^2
=d_dt(dtheta1_dt)
=d_dt(theta2)
d^2theta_dt^2=dtheta2_dt...........(3)
So, dtheta2_dt+(b/m)*theta2+(g/l)*sin(theta1)=0....(4)
dtheta1_dt=theta2........(5)
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1).......(6)
Now we have two equations,lets club into matrix form.
d_dt[theta1; theta2]=[theta2; -(b/m)*theta2-(g/L)*sin(theta1)]
and dtheta_dt=[dtheta1_dt; dtheta2_dt]
To find dtheta_dt matrix, we will create ODE function
Note: theta is the solution array
Please find below ODE function as name given "ode_function"
function [dtheta_dt]=ode_function(b,g,l,m,t,theta) 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
In above this function we have to find values for result. Result has two columns, 1st for displacement and 2nd velocity. Thats why i am assigning,
theta1=theta(1) and theta2=theta(2)
The result part we will cover next,
Now, we will find values for displacement and velocity w.r.t. time in a matrix form of table result.
Then call above ode_function into main program,
close all clear all clc
b=0.05; l=1; m=1; g=9.81;
theta_0=[0;4]; t_span=linspace(0,20,300); [t,result]=ode45(@(t,theta) ode_function(b,g,l,m,t,theta),t_span,theta_0);
figure(1) hold on plot(t,result(:,1),'color','r'); plot(t,result(:,2),'color','b');
xlabel('Time')
ylabel('Displacemt/Velocity')
legend('Displacement','Velocity')
grid on
In ode45 synatx we defining value of b,g,l,m are assumed and then we define time span t_span and initial conditions by theta_0
The result of above progrm as belows:
In order to create animation vbideo i used counter (ct) command and movie command.
M is the variable to store all the output result figures and use for loop to complete time span length(t) i.e. 300 intervals.
close all clear all clc
b=0.05; l=1; m=1; g=9.81;
theta_0=[0;4]; t_span=linspace(0,20,300); [t,result]=ode45(@(t,theta) ode_function(b,g,l,m,t,theta),t_span,theta_0);
figure(1) hold on plot(t,result(:,1),'color','r'); plot(t,result(:,2),'color','b');
xlabel('Time')
ylabel('Displacemt/Velocity')
legend('Displacement','Velocity')
grid on
ct=1; for i=1:length(t)
x0=0; y0=0;
x1=l*sin(result(i,1));
y1=-l*cos(result(i,1));
figure(2)
plot([-2 2],[0 0],'linewidth',3,'color','r');
line([x0 x1],[y0 y1],'linewidth',4);
axis([-2 2 -2 2])
hold on
plot(x1,y1,'o','markers',20,'markerfacecolor','g');
hold off
grid on
pause(0.03);
M(ct)=getframe(gcf)
ct=ct+1;
end
movie(M)
videofile=VideoWriter('Second_order_ODE.avi','uncompressed avi')
open(videofile)
writeVideo(videofile,M)
close(videofile)
Please find below link for Animation video of Simple Pendulum Single ODE
https://youtu.be/-mjA1DrRmgY