Solving second order ODEs using MATLAB

In this project we will be solving ODE equation with MATLAB, The ODE is:

Inputs for this project are:

L=1 metre,

m=1 kg,

b=0.05,

g=9.81 m/s2.

angular displacement=0,angular velocity=3 rad/sec at time t=0 i.e. theta_0=[0;3]

time span is 0 to 20. t_span=linspace(0,20,200). I have used 200 intervals.

Now to solve ODE, we first simplify it:

Assume theta=theta1,

Therefore, dtheta_dt=dtheta1_dt

Now, dtheta1_dt=theta2 --------------------- (1)

d^2theta1_dt^2=dtheta_dt

dtheta2_dt+(b/m)*theta2+(g/L)*sin(theta1)=0

dtheta2_dt=-(b/m)*theta2-(g/L)*sin(theta1) ---------------- (2)

Now we have two equations, which we will put intomatrix form,

d_dt[theta1; theta2]=[theta2; -(b/m)*theta2-(g/L)*sin(theta1)]

and

dtheta_dt=[dtheta1_dt; dtheta2_dt]

Now to find dtheta_dt matrix, we will create a function, because we need to find for 200 times as we have taken 200 values for t.

Check the code for function,

function dtheta_dt=ode_function(b,g,L,m,theta,t)
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

here I have used ode_function as name for the function:

Explaination:

• In this function we have to find values for result, the result part I will cover in next code. Result has two columns, 1st for displacement and 2nd velocity. Thats why we have assigned theta1=theta(1) and theta2=theta(2).
• We have provided equations to find dtheta1_dt and dtheta2_dt.b,g,L,m,t,theta are the variables used in the function so they are put in the bracket.

Now, we will find values for displacement and velocity w.r.t. time in a matrix form of table result.

We will use ode45 function to solve ODE in MATLAB.

Check code:

close all 
clear all
clc

%Inputs 
L=1 ;
m=1 ;
b=0.05;
g=9.81;

%Initial conditions, at t=0,
theta_0=[0;3];

%time interval is of20 sec. 
t_span=linspace(0,20,200);

%solving ode
[t,result]=ode45(@(t,theta) ode_function(b,g,L,m,theta,t), t_span, theta_0);

figure(1);
plot(t,result(:,1),'color','b')
hold on 
plot(t,result(:,2),'color','r')
xlabel('time')
ylabel('displacement/velocity')
legend('displacement', 'velocity')
grid on

You can see, I have used holdon command to put the results of displacement and velocity w.r.t. to time in one graph. I have named this as figure(1).

Explaination:

• MATLAB uses ode45 function to solve ODE.
• In the ode45, first we define function handles using @(). In this, ode function handles aret and theta.ode_function is used in ode45 solver to solve 200 intervals. t_span is time interval which is 0 to 20 sec. theta_0 is inital condition at t=0.

Check the graph:

I have used legend command to differentiate displacement and velocity.

Now, We have to make an animation of this, For this, I have used counter command and movie command. I have used M variable to store all the figures. For loop is used to generate figures with respect to the lenght(t).

Check the final command:

close all 
clear all
clc

%Inputs 
L=1 ;
m=1 ;
b=0.05;
g=9.81;

%Initial conditions, at t=0,
theta_0=[0;3];

%time interval is of20 sec. 
t_span=linspace(0,20,200);

%solving ode
[t,result]=ode45(@(t,theta) ode_function(b,g,L,m,theta,t), t_span, theta_0);

figure(1);
plot(t,result(:,1),'color','b')
hold on 
plot(t,result(:,2),'color','r')
xlabel('time')
ylabel('displacement/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',2)
    
    line([x0 x1],[y0 y1],'linewidth',4)
    hold on 
    plot(x1,y1,'o','markers',20,'markerfacecolor','r')
    hold off 
    axis([-2 2 -2 2])
    grid on 
    pause(0.03)
    M(ct)=getframe(gcf)
    ct=ct+1
end

movie(M)
videofile=VideoWriter('week_three_challenge','uncompressed avi')
open(videofile)
writeVideo(videofile,M)
close(videofile)

I have set the axis so that in animation the size of the frame remain unchanged.

Error Faced:

1. I have used theta_0 before t_span, which does not fit into ode45 solver's code. Time span must come before intial condition. Error that I got was the size of the matrix is not equal.
2. I have used cosine to find x1 and sine to fin y1 which lead motion of the pendulum wrong. Later I figured it out. Please check code.
3. I have wrongly used VideoWriter command. I had used videowriter which is a typographic error.

I have uploaded animation of the ODE, please check link:

Done.