Week 3 - Solving second order ODEs

AIM: -

To solve the ODE using MATLAB and plot the graph.

GOVERNING EQUATIONS: -

The equation below is a 2^nd order differential equation that represents the damped motion of a pendulum.

`(d^2 θ)/(dt^2 )+ b/m dθ/dt+ g/L sinθ=0`

Where,

`theta` = angular displacement (radian)

b = damping coefficient

m = mass of the suspended body (kg)

L = Length of the suspended bob string (m)

t = time (sec)

g = acceleration due to gravity (m²/s)

OBJECTIVE: -

Followings are the objective of the project:

• Solve the ODE using MATLAB.
• Plot the displacement and Velocity graph using MATLAB .
• Create an animation of Pendulum motion for the specific period of time.

BODY: -

Following is the process which is used by MATLAB to solve ODEs:

Let `θ= θ_1` And

`(dθ_1)/dt= θ_2`

so the above differential equation becomes –

`(dθ_2)/dt+ (b/m)θ_2+ ( g/L)sinθ_1=0`

• In the Matrix form the above equations can be written as

`d/dt [[θ_1],[θ_2 ]]=[[θ_2],[- (b/m)θ_2-(g/L)sinθ_1 ]]`

# CODE – 1

Firstly, the constants are defined and then the value of some constant are provided like L, m, b, g. then we have defined theta_0 is a 2-by-1 matrix where the first row define the displacement and the second row define angular velocity . t_span shows the time period of solution of the solution which is 20 sec and we have divided it into 100 equal divisions.

The ODE solution format is as given below –

[independent_variable, result] = ode45(@ (independent_variable, dependent_variable)

ode_function_to_call, solution_span, initial_value_matrix);

Here,

• Independent_variable is time (t).
• Result is matrix with two columns 1^st one for Displacement and 2^nd one is for velocity.
• ode45 is the solver for solving ordinary differential equation in MATLAB.
• @ (independent_var, dependent_var) is used for defining the variable of the equation.
• ode_function_to_call is used to call the ode function which is in the format of equation (2).
• solution_span determines the span to calculate the solution as specified earlier.
• initial_value_matrix the values contains the initial values used to solve the equation.

In Next step we have use for loop used the generated frame for animation. For this the frame is initially initialised and after that the loop is started. Here, the for loop is used with i = 1 to length of 1^st column of matrix result we have got using ode45.

The rotational axis of the pendulum is fixed in each loop to (0,0) and the co-ordinate of the bob is obtained as (-L*sin(result(i,1)), -L*cos(result(i,1)).Here we have used negative sign as we know lower axes are by default negative so to show movement of bob below the origin we have used negative sign. Here i takes the value of 1^st column of the result we have got by solving ode.

Next we have use subplot command, which divide the plot area into several sections to plot animation as well as the variation of displacement and velocity of bob w.r.t time in same frame which would help us to understand the variation of velocity with the location of bob.

subplot(6,2,[1:4]) divides the plot section into 6 rows and 2 column i.e. total of 12 small subplots of which subplots 1-to-4 is combined into a single plot area in which we have plot the movement of pendulum. The command below of it plots the animation of the pendulum. plot([x0,x1],[y0,y1],'linewidth',1,'color','k') plots the suspending massless thread of width 1 unit and of black colour. Using hold on to hold the subplot and then using plot(x1,y1,'-o','markers',10,'markerfacecolor','k') plots the bob with black colour and of size 10 units and then we have fixed the axis to have a better view of the animation.

The next subplot subplot(6,2,[7 8]) the displacement graph is plotted and with plot(t(i),result(i,1),'o','markers',5,'markerfacecolor','g') the small marker is overlapped on the displacement graph using hold on command that shows the real time position of on overall curve. Here we have used 7^th and 8^th block instead of 5^th and 6^th as we wanted to leave space for titles otherwise it would have overlapped with other graph.

The next subplot subplot(6,2,[11 12]) is the same as displacement plot but with Velocity as variable for which we have selected second column of the solution we have got from results of ode45.

The pause() command pauses the loop for the specified time interval and finally using getframe(gcf) the frame is captured in animation vector. Finally, by increasing the frame with the next one the loop ends.

A videofile named PendulumAnimaion.avi is created using VideoWrite command. It is then opend using open command, using writeVideo() command the frames stored in animation vector is written to the videofile and finally the videofile is closed using close() command.

# CODE – 2

To simplify our code we have defined our differential equation at different location and called it in 1^st code.

The code-2 is the ODE Function code (ode_func) which is called by the ODE solver (ode45). The input parameters are – t, theta, L, m, b, g.

The theta matrix on the LHS has the values theta1 and theta2 equal to 1^st and 2^nd element of theta matrix i.e. theta(1), theta(2).

The 1^st equation is represented by d_theta1_dt = theta2, and the other one by d_theta2_dt = -((b/m)*theta2) - (g/L)*sin(theta1).

Finally the equation ends with the output variable d_theta_dt as a column vector with two rows having the above equations respectively.

CONCLUSION: -

The above code can used as an example to solve the ODEs of any order. The equations in the ode calling function should be in correct order. The final output variable of the ode function must be a column vector.

There might be some confusion in deciding the axis of the of the animation plot which can be corrected by hit and trail method i.e. running the code each time and then deciding the best fit axis.

#LINK

https://youtu.be/QxnMdUHUxmk

CODE-1

close all
clear all
clc

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

theta_0 = [0;3];

t_span = linspace(0,20,100);

[t, result] = ode45(@(t,theta) ode_func(t,theta,L,m,b,g),t_span,theta_0);

frame = 1;


for i = 1:length(result(:,1))
    
    x0 = 0;
    y0 = 0;
    x1 = -L*sin(result(i,1));
    y1 = -L*cos(result(i,1));
    
    subplot(6,2,[1:4])
    plot([x0,x1],[y0,y1],'linewidth',1,'color','k')
    hold on
    plot(x1,y1,'-o','markers',10,'markerfacecolor','k')
    axis ([-1,1,-1.5,1.5]);
    title('PENDULUM')
    hold off
    
    subplot(6,2,[7 8])
    plot(t,result(:,1))
    hold on
    plot(t(i),result(i,1),'o','markers',5,'markerfacecolor','g')
    hold off
    title('DISPLACEMENT')
    
    subplot(6,2,[11 12])
    plot(t,result(:,2))
    hold on
    plot(t(i),result(i,2),'o','markers',5,'markerfacecolor','g')
    hold off
    title('VELOCITY')
    
    pause(0.01)
    
    animation(frame) = getframe(gcf);
    frame = frame + 1;
end 

videofile = VideoWriter('PendulumAnimation.avi','Uncompressed AVI');
open(videofile);
writeVideo(videofile,animation)
close(videofile)

CODE-2
function [d_theta_dt] = ode_func(t,Theta,L,m,b,g)

Theta1 = Theta(1);
Theta2 = Theta(2);
d_theta1_dt = Theta2;
d_theta2_dt = -((b/m)*Theta2) - (g/L)*sin(Theta1);
d_theta_dt = [d_theta1_dt; d_theta2_dt];

end