Week 3 - Solving second order ODEs
MATLAB program to simulate the motion of a simple pendulum.
Aim: To plot a graph of velocity and displacement versus time using ode45 solver, and hence to create an animation of the simulation of a simple pendulum.
Theory:In mathematics, an ordinary differential equation (ODE) is a differential equation containing one or more functions of one independent variable and the derivatives of those functions.A Second order differential equation is an equation involving the unknown function y, its derivatives y' and y'', and the variable x.
Governing Equation:
`(d^2 theta )/dt^2+b/m*(d theta)/dt+g/l*sin(theta)=0`
where
- theta - the angular dispacement (in degree)
- b - damping coefficient
- m - mass of the bob (in kgs)
- l - length of the string (in m)
- g - acceleration due to gravity (in `m/s^2`)
Matlab code:
1.Main code
%program to simulate the movement of a simple pendulum
clc
clear all
close all
b=0.05; %Damping coefficient
l=1; %length of the string
m=1; %mass of the bob
g=9.81; %acceleration due to gravity
%initial conditions
theta_0=[0;3]; %determines the initial displacement and velocity
t_span=linspace(0,20,500); %time span of 20 seconds
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,l,m,g),t_span,theta_0); %defining ode45
figure(1) %first figure
hold on
plot(t,results(:,1)) %plot of displacement vs time
plot(t,results(:,2)) %plot of velocity vs time
axis([0 20 -5 5]); %axis specification
xlabel('time') %labelling the x axis
ylabel('Angular displacement and velocity') %labelling the y axis
legend('Angular Displacements','Angular velocity') %indicates the graph
title('Angular velocity and displacements v/s time graph') %title of the graph
hold on
grid on
%animation
figure(2)
ct=1; %counter variable
for i=1:length(results(:,1)) %begining of the for loop(animation)
theta=results(i,1); %value of theta of i^th row and first column
%determination of the coordinates
x0=0;
y0=0;
x1=-3*sin(theta);
y1=-3*cos(theta);
%plotting the results
plot([-1,1],[0,0],'linewidth',3,'color','g') %plotting the slab
hold on
plot([x0,x1],[y0,y1],'linewidth',4,'color','b') %plotting the string
hold on
plot(x1,y1,'O','markers',30,'markerfacecolor','r') %plotting the bob using markers command
axis([-5 5 -5 5]); %axis specification
pause(0.01)
title('simulation of pendulum') %title of the simulation
M(ct)=getframe(gcf); %collecting each frame using getframe command
hold off
ct=ct+1; %incrementing the counter variable by one
end %end of the fo loop
movie(M) %creating a movie
videofile=VideoWriter('ode_code.avi','Motion JPEG AVI'); %videofile is created
open(videofile) %opening the videofile
writeVideo(videofile,M) %writing all the frames to videofile
close(videofile) %closing the videofile
Code Explanation:
- The code is started with clc command which clears the command window before initiating the program.
- The values of the constants in the differential equation are determined
L=1 metre,
m=1 kg,
b=0.05.
g=9.81 m/s2.
- The initial conditions such as initial displacement and initial velocity is specified as column vector named [theta_0].
- The span of time is determined using linspace command which generates 500 values between 0 to 20.
- The ode45 solver is used inorder to solve the second order differential equation.
- The displacements and the velocities for a span of 20 seconds are stored in the form of array inside results.
- Using plot command the graph of displcement and velocity versus time is drawn.
- The xlabel and ylabel is specified
- The title of the graph is also determined using title command
- Inorder to distinguish between the two graphs legend command is used. The application of this command is reflected at the top right corner of the graph.
- All these plots are drawn in figure 1.
- The figure 2 is mainly used for the animation pupose.
- The counter variable is specified and the for loop begins for the purpose of the animation.
- The value of i is specified using the length of the displacement values for a span of 20 seconds.
- The coordinates of the string is also specified.
- The bob is indicated using the markers command and its color is selected using markersfacecolor command.
- The title of each plot is kept same so that the same gets displayed in the animation.
- Each plot is considered as a frame and are dumped into an array named M.
- A movie is created using all the frames.
- Finally an animation is created using VideoWriter command which stiches these frames together.
2. Function code
%function of pendulum simulation
function [dtheta_dt]=ode_func(t,theta,b,l,m,g) %defining the function
theta1=theta(1); %theta1 is the first value of theta
theta2=theta(2); %theta2 is the second value of theta
dtheta1_dt=theta2; %first derivative
dtheta2_dt=-(b/m)*theta2-(g/l)*sin(theta1); % second derivative
dtheta_dt=[dtheta1_dt;dtheta2_dt]; array of the two derivatives
endExplanation of the code:
- The function is saved as ode_func
- an array of dtheta/dt has to be formed which contains theta1 and theta2 as the values of theta
- theta1 will be the first value of theta
- theta2 is the second value of theta
- dtheta_dt is a function of t,theta,b,l,m,g
- The values of dtheta1_dt and dtheta2_dt are dtermined using above formulas
- These two values will be stored in the array [dtheta_dt]
- The same values are taken while plotting the graph.
Explanation of ode45:
ode45 is a solver which converts the higher order differential equations into a differential equations of the order one, and solves it.
MATLAB’s standard solver for ordinary differential equations (ODEs) is the function
ode45. This function implements a Runge-Kutta method for efficient computation.It is designed to handle the following general problem:
dx/dt = f(t, x), x(t0) = x0, where t0=0.
where t is the independent variable,
x is a vector of dependent variables to be found and
f(t, x) is a function of t and x.
The mathematical problem is specified when the vector
of functions on the right-hand side f(t, x), is set and the initial conditions,
x = x0 at time t0 are given.
Syntax of ode45:
Where
- odefun is the name of the function . This is the function where we will input the system of first
order ode’s to be integrated. - tspan is the vector defining the beginning and end limits of integration, as well
as how large we want our time steps to be. i.e. if we are integrating from
t = 0 to t = 20, and want to take 500 time steps, then
tspan=linspace(0,20,500) - y0 is the initial condition .
The ode45 code used in the program:
[t,results]=ode45(@(t,theta) ode_func(t,theta,b,l,m,g),t_span,theta_0);
- [t,results] - This is array which contains values of time, displacements and velocity
- ode45 - calling the solver ode45
- (t,theta) -It tells that it cotains dependent variable theta and independent variable t.
- ode_func(t,theta,b,l,m,g)-The function ode_func is called which contains all the input values specified in the bracket.
- t_span- initial condition of time is spcified.
- theta_0- This is the initial displacement and velocity.
Output
There are two outputs in this program
- figure(1)-The graph of Angular dispacements and Angular velocities versus time
- The animation of the simulation of the pendulum.
Angular velocity and Angular displacement versus time
The link of the animation is provided below:
https://drive.google.com/file/d/1uc8mItbheGcFAQNk_VWxQcPjwf46a1ur/view?usp=sharing
Errors faced:
Basic typing errors:
Error shown in command window
Typing error in line 52: 'linewith' was the error , which was resolved later.
Conclusion:
ode45 is a very good parameter that can be used to solve the higher order differential equation in MATLAB.
Reference: