Matlab Program to solve Ordinary Differential Equation of 2nd order using ODE45 Matlab Solver

• Aim: To create a program for solving a second order ordinary differential equation using ODE45 solver.
• Equation and Expressions: The following equations are used for calculating and one as an example:

1. `(d^2(theta))/dt^2+b/m*(d(theta))/dt+g/l*sintheta=0`

This is the example we are taking for solving, using the Matlab solver ODE45. Similarly we can change the input parameters for a different expression for evaluating the results. The above expression is for motion of a simple pendulum.

Where `theta` is the angular displacement

b is the damping cofficient

m is the mass of the pendulum in Kg

l is the lenght of the pendulum in meter(m)

g is the acceleration of gravity in m/`s^2`

2. `d^2=(x_2-x_1)^2+(y_2-y_1)^2`

This relation from the pythagoras theorem is used to calculate the y-coordinate of the pendulum as it follows the oscillating path whereas the x-coordinate is calculated from solving the differential equation.

Where d is the distance between the two coordinates

`x_1` and `y_1` are the coordinates of the fixed point of the pendulum typically for this example(0,1)

`x_2` and `y_2` are the coordinates of the pendulum itself

• Objective: We create a Matlab program to solve an 2nd order ODE using the ODE45 solver to analyse the results. For the above example over a period of time the `theta` is going to vary according to the differential equation which then will provide us the angular displacement and velocity for that specified time. Similarly we can analyse the behaviour of any other ordinary differential equation.

Workflow: Firstly we define the parameters of the pendulum which are b,l,m then we need to define the limit of time over which the output is need to be taken with the initial conditions of the differential equation.

Now we use the ODE45 solver to get the solution of the differential equation, from there we use the output data to plot the angular displacement of the pendulum.

So to make an animation of the different angular position of the pendulum over the time which is our output also, we use a loop for plotting each output position and capturing each frame. After that we use the movie command to stitch the all these captured frame.

• Matlab Program:

The language used is Matlab Scripting language.

%Program to solve 2nd order differential equation using ODE45 solver
clear global;
clc;

%Inputs parameters of the pendulum
l=1;     %lenght in meters
m=1;     %mass in kg
g=9.81;  %acceleration of gravity in m/s^2
b=0.3;   %damping cofficient

%Inputs for the Ode funtion
angle=[0;3];              %initial conditions at time t=0;
t=linspace(0,20,200);     %time for which integeration will takes place

%solving the differential expression using the ode45 solver
[ty, output_x]=ode45(@(t,angle) ode_func(t,angle,b,g,l,m),t,angle);

%calculating the variation in vertical while pendulum oscilates
%horizontally, using d^2=(x2-x1)^2+(y2-y1)^2
output_y=1-(l^2.-(output_x(:,1)).^2).^(1/2);

% plotting the outputs
for i=1:200
plot([output_x(i,1) 0],[output_y(i,1) 1],'-o','linewidth',2,'MarkerSize',15,'MarkerEdgeColor',[1 0 0.1],'MarkerFaceColor',[1 0.5 0]);

hold on;
plot([0 0],[-1.5 1],'-.','MarkerIndices',1:50:5);
hold off;
%plot(output(i,2),1,'-o');
xlabel('positon');
axis([-1.5 1.5 -0.5 1.5]);
%capturing each frame using getframe and then putting in it variable fig
fig(i)=getframe(gcf);
end
    movie(fig);
    %Creating a avi file providing its name and profile 
    video = VideoWriter('Kinematics_of_simple_pendulum.avi','Uncompressed AVI');
    open(video);
    %writing the frames to the avi file
    writeVideo(video,fig);
    close(video);

Explanation:

1. We take the differential equation which defines the motion of a simple pendulum, a 2nd order differential equation. So for a simple pendulum there we take the inputs of its lenght, mass, damping cofficient and the acceleration of gravity(constant).
2. Then we need to define the time span over which we want the results. These are basically the limits of the intergration over which the solution of the differential equation will be taken. We must provide the time because the solver will provide the variation in the angular dispalcement and velocity with variation in time. So over what time from the starting and ending need to be stated. We need to provide the initial conditions of the differential equation for the starting time, with intial angular displacement be 0 and velocity be 3.
3. ODE45: It is a 1st order ordinary differential equation solver which is based on the 4th order Runge Kutta method. Unlike other solvers in it we need not to define the step size because it has an adaptive step size to keep the error in tolernace level. We need to provide the time span and also the initial conditions for the equation at time t=0 as it is based on R. K. method.
4. Syntax: There will be two outputs one is the time with respect to which the second output angular displacement and velocity is calculated. The output_x is a column vector having two columns, the first one is angular displacement and the second one is velocity. We use the ode45 command and describe an anonymous function having first the independent variable then the dependent variable declared which in this case are the time and angle. Further this we need to define the differential equation for which we create a function and call that here namely(ode_func). Now according to the need for solver we need to provide the time span over which integration will be done and then the initial conditions, careful the order needs to be same as time then conditions.
5. The output_x is the variation of the pendulum in horizontal direction with respect to the fixed point(0,1), in which we are using the displacement column vector. So to find out the variation of the pendulum in the vertical direction we use the expression from the pythagoras theorem as the lenght of the pendulum remains constant from the fixed point to the pendulum at each position while it oscillates. So it will also be a column vector(output_y) containing the values for the variation of the pendulum in vertical(y) direction with respect to each value for pendulum in horizontal(x) direction.
6. Now comes the plotting of the coordinates of the pedulum, to plot each position of the pendulum one after another we need to use a FOR loop so that we can plot the coordinates sequentially. We use the MarkerIndices command to create circle at the coordinates of the pendulum.Then to color the faces and edges the marker we use the MarkerEdgeColor and MarkerFaceColor.
7. Now to create a marked line passing through the centre we use another plot command but to keep both the plot in the same figure window we use hold on commad. After which we need to hold off command also otherwise after each iteration the previous position of the pendulum will also be there.
8. For creating an animation of the pendulum motion we capture each frame carrying each position using the getframe command and store these frames in a variable 'fig'.We use the function 'movie' to create an animation of the captured frames.
9. We need to create the file name and type of file for the animation so we use the command VideoWriter to create an AVI file. Now we open this file write the frames using the command writeVideo and then closes it.

Function ode_func:

function [dangle]=ode_func(~,theta,b,g,l,m)

theta1=theta(1);
theta2=theta(2);

dtheta1=theta2;
dtheta2=-(b/m)*theta2-(g/l)*sin(theta1);

[dangle]=[dtheta1; dtheta2];
end

Explanation:

We create this function in which we describe the differential equation which we need to solve.But we are using the ODE45 solver which solves/takes only 1st order differential equation, so we use a method to break it down into two 1st order differential equations. We assume the change in dependent variable as another variable providing us the first one while putting this newly assumed variable in the original equation will give us the second one, much similar to integration by substitution. So assuming `theta` to be `theta_1` and then `(d(theta_1))/dt=theta_2` we get the required equations. The output of the function needs to be a column vector.

• Ouput:
  https://drive.google.com/file/d/1te2IARbNeJ7M-eV1Lnkk52-M70hpxZFm/view?usp=sharing

https://drive.google.com/file/d/1S3gyTX2ZgQw7JDQhdyH6C-5kAD7mkkUz/view?usp=sharing