Week 3 - Solving second order ODEs
MAIN SCRIPT
% Program to solve the second order ODE using build in function ode45 and to create an animation using the solution of the equation
% 2nd-order equation: d^2(theta)/dt^2 + (b/m)*d(theta)/dt + (g/L)*sin(theta) = 0
% substituting theta = theta1
% substituting d(theta)/dt = d(theta1)/dt = theta2
% by making these substitutions we can rewrite the 2nd-order ODE as two first-order ODEs
% d(theta1)/dt = theta2 ------------ 1st eqn
% d(theta2)/dt = -(b/m)*(theta2)-(g/L)*sin(theta1) ---------- 2nd eqn
clear
close all
clc
% Initializing the variables for input values
L = 1; % Length of pendulum (m)
m = 1; % mass of the pendulum (kg)
b = 0.05; % damping coefficient (dimensionless)
g = 9.81; % acceleration due to gravity (m/s^2)
ct = 1; % initializing the counter
% Specifying the time interval and initializing the variable t_span for it.
t_span = linspace(0,20,500); % time interval - 0 to 20 sec, with 500 data points
% Specifying the initial conditions at time t = 0 sec.
% initial angular displacement = 0
% initial angular velocity = 3 rad/s
theta_0 = [0;3]; % column vector representing initial conditions at t = 0.
% calling ode45 function to solve the ODE. Two outputs are returned by the function
[t,solutions] = ode45(@(t,theta) func_ode(t,theta,L,m,b,g),t_span,theta_0);
% 't' is the time column vector
% 'Solutions' is an array with two columns.
% The first column of array solutions represents angular displacement with respect to time
% The second column of array solutions represents angular velocity with respect to time
% plotting
figure(1) % 1st figure will have plots for time,angular displacement and angular velocity
plot(t,solutions(:,1)) % time and angular displacement plot.
hold on % command hold on holds the plot to the same figure
plot(t,solutions(:,2)) % time and angular velocity plot.
xlabel('Time') % labelling horizontal axis as time
ylabel('Solutions') % labelling vertical axis as Solutions
title('Solution of Simple Pendulum Equation with damper') % giving title to the plot
legend('Angular Displacement (rad)','Angular Velocity (rad/s)') % creating legend to show the lines used for displacement and velocity
% Creating the animation movie for pendulum
ang_disp = solutions(:,1); % assigning values of angular dispacement from the array solution to a variable ang_disp. This will create a column vector ang_disp.
% using for loop to get the series of plots to create an animation
for i = 1:length(ang_disp)
% Specifying the coordinates of the fixed end of the pendulum
xf = 0;
yf = 0;
% calculating the coordinates of the end with bob.
xb = L*sin(ang_disp(i));
yb = -L*cos(ang_disp(i));
figure(2) % new figure for the animation plot
plot([-1 1],[0 0],'color','k','linewidth',6) % plotting the line representing fixed surface
hold on
plot([0 0],[-1.5 1.5],'-.','color','k','linewidth',1) % plotting the line representing mean position
hold on
plot([xf xb],[yf yb],'color','k','linewidth',4) % plotting the line representing pendulum rod
hold on
plot(xf,yf,'o','markersize',5,'markerfacecolor','w') % plotting the cirle marker representing the fixed end joint between pendulum and fixed surface
hold on
plot(xb,yb,'o','markersize',20,'markerfacecolor','c','markeredgecolor','k') % plotting the circle marker representing the bob
hold off % hold off so that the next frame does not carry this plot.
axis([-1.5 1.5 -2 0.5]) % scaling the x and y axes
xlabel('X-axis') % labelling x axis
ylabel('Y-axis') % labelling y axis
M(ct) = getframe(gcf); % funciton getframe returns a movie frame. It takes snapshots and store them in the variable M
ct = ct+1; % counter variable ct will help changing the indices of the M array with each loop
end
movie(M) % Playing the movie
% creating and saving the video file
videofile = VideoWriter('Challenge_Simple_Pendulum.avi','Uncompressed AVI'); % creating a file in which the movie can be written.
open(videofile) % openning the file
writeVideo(videofile,M) % writing data to the file
close(videofile) % closing the file
User-defined function - func_ode
function [dtheta_dt] = func_ode(t,theta,L,m,b,g)
% this function defines the system of equations which is created after rewriting the second-order ODE as a first-order equation system
% d(theta1)/dt = theta2 ------------ 1st eqn
% d(theta2)/dt = -(b/m)*(theta2) - (g/L)*sin(theta1) ---------- 2nd eqn
% initializing variables to assign initial condition values
theta1 = theta(1); % angular displacement at (t = 0) = 0
theta2 = theta(2); % angular velocity at (t = 0) = 3 rad/s
% initializing the variables and assign the first order equations to them
dtheta1_dt = theta(2); % 1st eqn
dtheta2_dt = -(b/m)*theta2 - (g/L)*sin(theta1); % 2nd eqn
% initializing the output array (dtheta_dt) and should return 2x1 column vector.
[dtheta_dt] = [dtheta1_dt; dtheta2_dt];
endThe objective of this program is:
1) To solve the second order ODE using built in function ode45.
2) To create an animation using the solution of the equation.
This 2nd-order ODE represents the motion of a simple pendulum with damper:
d^2(theta)/dt^2 + (b/m)*d(theta)/dt + (g/L)*sin(theta) = 0
substituting theta = theta1
substituting d(theta)/dt = d(theta1)/dt = theta2
by making these substitutions we can rewrite the 2nd-order ODE as two first-order ODEs
d(theta1)/dt = theta2 ------------ 1st eqn
d(theta2)/dt = -(b/m)*(theta2)-(g/L)*sin(theta1) ---------- 2nd eqn
Explanation of the main script
First I've initialized the variable for input values.
L is the length of pendulum in meter, m is the mass in kg, b is the damping coefficient (dimensionless) and g is acceleration due to gravity (m/s^2)
L = 1;
m = 1;
b = 0.05;
g = 9.81;
Then I've initialized a counter variable which will be required later.
ct = 1;
Then I've used linspace command to specify the time interval from 0 to 20 seconds with 500 data points.
t_span = linspace(0,20,500);
The given initial conditions at time (t=0) are:
Initial angular displacement = 0
Initial angular velocity = 3 rad/s
For initial conditions I've used a column vector and assigned it to varable called theta_0.
theta_0 = [0;3];
Now to solve the equation I've used the build-in function (ode45)
To call the function I've used the syntax:
[t,solutions] = ode45(@(t,theta) func_ode(t,theta,L,m,b,g),t_span,theta_0);
Explaining the syntax
On the left side are the two output variables t and solutions. These variables will store the values returned by the ode45 function after solving the equation. On the right side the functions is called ode45 and inside the parenthesis @(t,theta) specifies the two main input variables. Then the user defined function that defines the given system of equation is specified and inside the parenthesis along with t and theta the rest of the input variables can be passed func_ode(t,theta,L,m,b,g). After this the time interval is specified t_span and lastly the initial conditions are specified theta_0
The function ode45 will solve the equations by integrating the system of equations for time interval t_span and initial conditions theta_0.
For each time value in vector 't' the function will return a corresponding value of angular displacement and angular velocity in the array 'solutions'
The output variable 't' is the time column vector
The output variable 'solutions' is an array with two columns. The first column represents angular displacement with respect to time and the second column represents angular velocity with respect to time.
After this I've plotted the time vs angular displacement and time vs angular velocity in figure 1. command hold on is used to plot both plots in same figure. I've labelled the x-axis as 'Time' and y-axis as volume. I've used title command to give the plot a title. To indicate the displacement and velocity plot I've used legend command.
figure(1)
plot(t,solutions(:,1))
hold on
plot(t,solutions(:,2))
xlabel('Time')
ylabel('Solutions')
title('Solution of Simple Pendulum Equation with damper')
legend('Angular Displacement (rad)','Angular Velocity (rad/s)')
To create the animation the coordinates of both the ends of pendulum are required. One end of the pendulum is fixed. The other end with the bob will move. The coordinates of the moving end with bob are calculated using the length of pendulum (L) and the values of angular displacements in the array solutions. I've assigned the values of angular displacement from the 'solutions' array to a variable called ang_disp.
ang_disp = solutions(:,1);
Now a for loop is used to calculate all the coordinates and plot them repeatedly to get frames for the animation. The loop variable i is initialized. It will go from 1 to the length of vector ang_disp which is 500 in this case. Meaning loop will run 500 times.
for i = 1:length(ang_disp)
Then I've specified the coordinates of the fixed end of the pendulum at (0,0) and assigned them to xf and yf
xf = 0;
yf = 0;
After this I've calculated the coordinates of the end with bob and assigned them to xb and yb. The x and y coordinated will be calculated for all the values of the angular displacements in the array ang_disp. The values of angular displacements will be provided by ang_disp(i) as the value of (i) changes from 1 to 500. Since the pendulum is fixed at 0,0 the moving end will lie on the -ve y-axis therefore the -ve sign is used for calculating yb.
xb = L*sin(ang_disp(i));
yb = -L*cos(ang_disp(i));
Now for the animation I've created a new figure
figure(2)
Now in this figure I've plotted all the plots required for the animation. The first plot will represent the fixed surface in black color and linewidth 6. The second plot represents the mean position shown by a dash dotted black line and linewidth 1. The third plot is the pendulum rod plotted in black and linewidth 4. The fourth plot is showing the joint between the fixed surface and the fixed end of the pendulum with a white circle marker having marker size 5. The fifth plot represents the bob in cyan color with a cyan circle marker having black markeredge and size of the marker is 20.
I've used hold on command for every plot so that all the plots are plotted in the same figure to shape a pendulum. After the fifth plot I've used hold off so that the coordinates do not over lap in the figure.
I've also specified the scales for x and y axes and also have labelled them
plot([-1 1],[0 0],'color','k','linewidth',6)
hold on
plot([0 0],[-1.5 1.5],'-.','color','k','linewidth',1)
hold on
plot([xf xb],[yf yb],'color','k','linewidth',4)
hold on
plot(xf,yf,'o','markersize',5,'markerfacecolor','w')
hold on
plot(xb,yb,'o','markersize',20,'markerfacecolor','c','markeredgecolor','k')
hold off
axis([-1.5 1.5 -2 0.5])
xlabel('X-axis')
ylabel('Y-axis')
The function getframe returns a movie frame and it is stored in the array M. The variable ct is used to specify the indices to store the frames. In this case there will be 500 frames stored inside the array M
M(ct) = getframe(gcf);
The counter ct is set inside the loop. Every time the loop runs the value of ct is incremented by 1.
ct = ct+1;
the for loop is closed here with an end command
end
Now to test I've played the movie using movie command
movie(M)
To save the movie I used the command VideoWriter and created an .avi file 'Challenge_Simple_Pendulum.avi', specified the format 'Uncompressed AVI' and assigned it to variable videofile. Then I opened the file using open command and wrote the data inside it using command writeVideo. Then I closed the file using close command
videofile = VideoWriter('Challenge_Simple_Pendulum.avi','Uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile).
Explanation of the user defined function - func_ode
The function func_ode defines the system of differential equations which is obtained by rewriting the 2nd-order ODE as 1st-order equation sytem.
d(theta1)/dt = theta2 ------------ 1st eqn
d(theta2)/dt = -(b/m)*(theta2)-(g/L)*sin(theta1) ---------- 2nd eqn
first I've initialized two variables theta1 and theta2 and assigned values of initial conditons to them.
theta1 = theta(1);
theta2 = theta(2);
Then I've initialized the variables for the equations
dtheta1_dt = theta(2);
dtheta2_dt = -(b/m)*theta2 - (g/L)*sin(theta1);
Then I've initialized the output array dtheta_dt. This array should be in the form of a 2x1 column vector.
[dtheta_dt] = [dtheta1_dt; dtheta2_dt];
Then I've closed the function
end.
This concludes the program.
There was one problem that I faced in the execution of the program. The use of hold on command was leading to the overlapping of plots with each iteration of loop and my animation was getting messed up. I fixed this by using hold off command. I've attached a screenshot of this problem. Apart from this I've attached the screenshot of workspace after execution of program to show all the variables created. Also, one screenshot of figure 1 and three screenshots of figure 2 are attached.
The animation video is uploaded to the youtube and here's the link: