Week 3 - Solving second order ODEs
ODE function program: pendulum animation
Matlab main code:
clear all
close all
clc
b=0.05;
g=9.81;
l=1;
m=1;
%initial condition
%given angular vel = 3 rad/s, linear velocity = r* angular velocity
theta_0 = [0;3];
%time points
t_span = linspace(0,20,500);
%solve Ode
[time,wave] = ode45(@(time,theta)functionode(theta,b,g,l,m), t_span, theta_0);
%assiging the displacement coloumn of waves to X
X=wave(:,1);
%opening a figure window in maximised state
figure1 = figure('WindowState','maximized');
ct = 1;
for i=1:length(X)
x2=X(i);
x1=0;
y1=0;
y2=((10)-(x2^2))^0.5;
y3 = -y2;
plot([x1 x2],[y1 y3],x2,y3,'.','markersize',50, 'linewidth', 2)
axis([-3 3 -3.5 1])
M(ct)= getframe(gcf);
ct = ct+1;
end
movie(M)
videofile = VideoWriter('pendulum','MPEG-4');
videofile.FrameRate = 25;
open(videofile)
writeVideo(videofile,M);
close(videofile)Matlab function code: functionode.m
function [dtheta_dt] = functionode(theta,b,g,l,m)
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
Output:
video link: pendulum video
Incase link dos not work, copy paste the url: https://drive.google.com/file/d/1LJ44-2cDOf5OYfvDU95-0wSn3d1z8Ogf/view?usp=sharing
Explanation:
Main code:
clear all
close all
clc
clearing previous data and command window for clean program exeution
b=0.05;
g=9.81;
l=1;
m=1;
Stating constant values of damping coefficient, gravitational acceleration, length of pendulum, and mass of pendulum bob.
%initial condition
%given angular vel = 3 rad/s, linear velocity = r* angular velocity
theta_0 = [0;3];
Theta_0 states the initial conditions required by the ode function. Displacement on left is zero and velocity on right is 3m/s. Angular velocity was given, which when multiplied with length or radius gives linear velocity.
%time points
t_span = linspace(0,20,500);
t_span is time duration and time steps for ode function. Linspace command divides the time duration of 20 seconds between 500 equal parts.
%solve Ode
[time,wave] = ode45(@(time,theta)functionode(theta,b,g,l,m), t_span, theta_0);
on the left hand side of the above equation, time and wave are the outputs of the ode function. Time here store the time data of values between 0 to 20 seconds, and wave stores both the displacement and velocity data in two columns.
The ode45 function is a generalised ode solver in matlab, and there are other ode solvers as well that are made for specific ode types, but ode45 is the most generalised on and can be used for any ode type.
Within the brackets of ode45 the function code ‘ functionode’ is called which is created by us, that involves the ode that we need to solve as shown in the function code. After the function, the time duration is either called or mentioned, and in this case it is called from a different location within the code. And the last variable theta_0 states the initial condition for the code.
%assiging the displacement coloumn of waves to X
X=wave(:,1);
when the ode is solved, the time and wave variables have the desired output of time, displacement and velocity values. But for creating an animation or plot of position of the pendulum, we desire only the displacement value and nothing else. Thus the variable ‘X’ is given all the values of the first column of the wave variable.
This way X is an array consisting the values of displacement.
%opening a figure window in maximised state
figure1 = figure('WindowState','maximized');
The default size of window of the figure is very small and matlab records screenshots of it, thus forming a very small sized video file. To have a better quality, this code opens up a figure window in the maximised state before any other function attempts to open it in default state.
ct = 1;
variable ct is used for counting the frames for forming the animation when the movie command calls on the array ‘m’ that stored all the frames of the plots.
for i=1:length(X)
For loop function is used to plot every single value of displacement, as it contains 500 values of it.
x2=X(i);
Here, x2 is stated as ith value of X array as the plot function cannot plot a complete array, an even if it does, it would beat the purpose of creating animation.
x1=0;
y1=0;
Stating the point where the pendulum is tied or hinged.
y2=((10)-(x2^2))^0.5;
Equation to find out the value of y co-ordinate for every x co-ordinate which represents the displacement from the centre position or origin. It is done with the help of equation of a circle in co-ordinate systems.
y3 = -y2;
here the value of y2 obtained is positive and thus above the origin of the graph and hence plots a reverse pendulum motion. To reverse the reversed plots, the y value is inverted to negative y and the plot appears to be normal.
plot([x1 x2],[y1 y3],x2,y3,'.','markersize',50, 'linewidth', 2)
This function plots the line as well as the solid marker that represents the bob.
Axis([-3 3 -3.5 1])
This function or code limits the axis of the graph so that the graph background is static and only the plotted line and points change giving the effect of and animation.
M(ct)= getframe(gcf);
This code collects the plotted graph screenshot and stores in the array ‘M’
ct = ct+1;
increasing the ct value by 1 so that next frame stored in a new space in ‘M’ and does not overwrite the previous frame.
End
For loop ends here.
movie(M)
videofile = VideoWriter('pendulum','MPEG-4');
videofile.FrameRate = 25;
open(videofile)
writeVideo(videofile,M);
close(videofile)
The above set of code initiates the process of creating the video file by fetching the frames stored in array ‘M’. The frame rate is set to 25 so that the actual speed of pendulum is obtained in the video.
Matlab function code: functionode.m
function [dtheta_dt] = functionode(theta,b,g,l,m)
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
The above code is a function code that has a set of equations that calculate the desired value. This function code can be called in any other program as long as they exist in the same folder. In this function code, the ODE mentioned in the queation is used to calculate the values.