Week 3 - Solving second order ODEs

AIM: -

To solve the ODE in MATLAB and also 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,

`θ`= angular displacement (radian)

b = damping coefficient (Ns/m)

m = mass of the suspended body (bob) (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:

• Write a program to solve the ODE in MATLAB.
• Plot the displacement and Velocity graph.
• Create an animation of Pendulum motion for the specified 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/Lsinθ_1=0` (1)

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

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

close all
clear all
clc

% Constants

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

% Initialisation

theta_0 = [0;3];

% Time span of result

t_span = linspace(0,20,500);

% ODE solution

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


% animation

frame = 1;

% Loop for graph and animation

for i = 1:length(result(:,1))
    
    x0 = 0;
    y0 = 0;
    x1 = -L*sin(result(i,1));
    y1 = -L*cos(result(i,1));
    
    subplot(3,2,[1:4])
    plot([x0,x1],[y0,y1],'linewidth',3,'color','b')
    hold on
    plot(x1,y1,'-o','markers',20,'markerfacecolor','g')
    axis ([-1,1,-1.5,0.5]);
    title('Pendulam Animation')
    hold off
    
    subplot(3,2,5)
    plot(t,result(:,1))
    hold on
    plot(t(i),result(i,1),'o','markers',5,'markerfacecolor','r')
    hold off
    title('Displacement')
    
    subplot(3,2,6)
    plot(t,result(:,2))
    hold on
    plot(t(i),result(i,2),'o','markers',5,'markerfacecolor','r')
    hold off
    title('Velocity')
    
    pause(0.4)
    
    animation(frame) = getframe(gcf);
    frame = frame + 1;
end 

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

The above code helps us to solve the above differential equation, plot the Displacement and Velocity graph and also the animation.

Firstly, the constants are defined and then the initial value required for solving is initialised. theta_0 is a 2-by-1 matrix with the rows as displacement and velocity initial values. t_span shows the time period of solution.

The ODE solution format is as given below –

[independent_var, result] = ode45(@ (independent_var, dependent_var) ode_function_to_call, solution_span, initial_value_matrix);

Here,

• Independent_variable is time
• Result is matrix 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.

Next step is to use for loop and capture the frame to generate the animation. For this the frame is initially initialised and then the loop is started. Here, the for loop is used with i = 1 to length of 1^st column of matrix result.

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))) where as specified earlier 1^st and 2^nd column of result matrix contains the displacement and velocity respectively at each timestep.

Using the subplot command, the plot is divided into several sections to plot animation as well as the variation of displacement and velocity of bob w.r.t time. subplot(3,2,[1:4]) divides the plot section into 3 rows and 2 column i.e. total of 6 small subplots of which subplots 1-to-4 is combined into a single plot area. The command below of it plots the animation of the pendulum. plot([x0,x1],[y0,y1],'linewidth',3,'color','b') plots the suspending massless thread. Using hold on to hold the subplot and then using plot(x1,y1,'-o','markers',20,'markerfacecolor','g') plots the bob with green colour and the fixing the axis to have a better view of the animation.

The next subplot subplot(3,2,5) the displacement graph is plotted and with plot(t(i),result(i,1),'o','markers',5,'markerfacecolor','r') the small marker is overlapped on the displacement graph using hold on command that shows the real time position of on overall curve.

The next subplot subplot(3,2,6) is the same as displacement plot but with Velocity as variable.

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 Pendulum.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.

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

The above code 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 matrix (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 `(dθ_1)/dt= θ_2` 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 is 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.

The link of the animation along with plot is provided below –

https://youtu.be/HMI0aWQkcyA