Solving Second order Ordinary Differential Equation

%solving 2nd order ode 
%by Roop Sannigrahi
close all
clear all
clc

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


%initial condition
theta_0 = [0;3];

%time points
t_span = linspace(0,20,200);

%solve ODE
[t,results] = ode45(@(t,theta) ode_func(t,theta,b,l,g,m),t_span,theta_0);

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('amplitude')
xlabel('time')
legend({'Angle (rad)','Angular speed (rad/s)'})

%% animation code



%for loop for animation
ct = 1;
for i=1:length(results(:,1))
    figure(2)
   
    %pendulum
    x0 = 0;    %giving fixed points
    y0 = 1;
    x1 = l*sin(results(i,1));
    y1 = -(l*cos(results(i,1)));
    plot([x0 x1],[y0 y1],'-o')
    axis([-1.5 1.5 -1.5 1.5]);     %Setting the limits of the plot
    pause(0.03)
    M(ct) = getframe(gcf);
     ct = ct +1 ;
    
end

movie(M)
videofile = VideoWriter('Pendulum_animation.avi','Uncompressed AVI');
open(videofile)
w riteVideo(videofile,M)
close(videofile)



function [dtheta_dt] = ode_func(t,theta,b,l,g,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

2^nd order ode solving program

In the program we are solving a 2^nd order ode used for the oscillatory motion of a simple pendulum.

The way the pendulum moves depends on the Newtons second law. When this law is written down, we get a second order Ordinary Differential Equation that describes the position of the "ball" w.r.t time.

In the above equation,

g = gravity in m/s2,

L = length of the pendulum in m,

m = mass of the ball in kg,

b=damping coefficient.

As this is a 2^nd order ode we will make it simpler by breaking it into two 1^st order ode namely

and .

So there are two 1^st order ode one for displacement and one for velocity.

So for the 1^st ode we need to provide initial conditions for displacement and for the second ode we need to provide initial conditions for velocity.

So we have taken theta_0 = [0;3] . here 0 is the starting position and 3 rad/sec is the starting angular velocity . We will integrate the ode equations to obtain a solution. So we have to define the value of time to integrate. So we use t_span and define it. 0-20 sec and 200 values in between using the linspace command.

Next there is a ode45 function , so it is an inbuilt MATLAB solver so the main inputs of this function is t and theta. There are other constants but these are the dependent and the independent variables.

Next we will define the function that we had made, so as said before we ill break the main 2^nd order ode into two 1^st order ode. The arguments that we will pass are time, theta, b , l, g ,m. We define them at variables namely dtheta1_dt and dtheta_2dt. As they are separated we need to combine them together and we make it in a matrix form.

Now we need to club it all together to get the result in [t,result] so this gives as 1 column matrix of time 200x1 and one 200x2 matrix the first column is of displacement and second column is for angular velocity. In the ode solver

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

Here we solve at time and the angular displacement theta as our ode is based on that primarily.

We are taking the plot of the solved ode , we will get the plot in a wave form as it is a oscillatory motion . we get the amplitude and time in the plot. We need to make the annimation of a pendulum so

Now we proceed to the animation , here we have used a simple for loop for animation and a counter variable ct that will help us later in the movie. So we take length of the angular displacement column for this as we will can the positon of the pendulum with theta that is angular displacement. We also need to use the velocity on the other column with respective displacement value so we use,

i=1:length(results(:,1)) this range. Here we get the above mentioned values. Now we need to the coordinates we have to convert it into cartesian form as they are in polar form . after that we did the plotting to get the animation as we are not using hold on we are getting an animation or I should say a cluster if images that is displayed with time.

For the movie we use the movie command and here the counter variable comes into play as this helps us to get each frame in M(ct) variable and then we increment ct after each loop.

Next we use the movie command to make the movie of the captured frames that we got using getframe() and we write the video file and then save it.

Output plot:

Animation link:

https://www.youtube.com/watch?v=Wg6p3ti1swk