SOLUTION OF A SECOND ORDER DIFFERENTIAL EQUATION USING LINEAR ALGEBRA AND VERIFYING IT WITH THE CALCULUS METHOD.

PROJECT TITLE: SOLUTION OF A SECOND ORDER DIFFERENTIAL EQUATION USING LINEAR ALGEBRA AND VERIFYING IT WITH THE CALCULUS METHOD.

OBJECTIVE: >Writing the code to solve the transient nature of pendulum and animating its motion.

>Simulating the motion for a definite time period starting from 0.

>Getting the position of the pendulum as a function of time and animating it.

The ODE which governs the pendulum motion is a second order differential equation. In this project, I have split the 2nd order differential equation into two coupled equations. Then both the equations are solved using numerical integration done by ode45 solver(4th order Runge Kutta method). The results for the displacement and the velocity are obtained and plotted against time.

To make it more interesting, I also used Linear Algebra for solving this system of Differential Equation. The concept of eigen values and eigen vectors is used to solve this equation and the results are validated with that of the numerical integration with good agreement.

CODE:

Here is the solution of the ODE using Matlab(Numerical Integration):

clear all
close all
clc

g = 9.81; % Acceleration due to gravity
l = 1;    % Length of Pendulum
m = 1;  % Mass of Pendulum
k=0.05    %damping coefficient   

% Initial Conditions, angular disp=0,angular velocity=3rad/s
theta_0 = [0;3];

% Time Points
t_span = linspace(0,20,600);


% Solve ODE
[t, results] = ode45(@(t,theta)ode_fun(t,theta,k,g,l,m),t_span,theta_0);

figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2),\'red\')
ylabel(\'Responses\')
xlabel(\'Time\')
legend(\'displacement vs time\',\'velocity vs time\')

figure(2)
ct=1;
for i=1:length(t_span)
       THETA = results(i,1);
    %Coordinates for pendulum
    x1 = l*sin(THETA);
    y1 = -l*cos(THETA);

     figure(2)
     % Plot string
     plot([0 x1],[0 y1],\'linewidth\',5)
     
     % Plot Support
     rectangle(\'Position\',[-0.5 0 1 0.5],\'EdgeColor\',\'k\',\'LineWidth\',3)
     
     %plotting of the ball
     ball=viscircles([x1,y1],0.2);
     grid on
     
     
      axis([-2 2 -2 1])
      title(\'Simple Harmonic Motion in Pendulum\')
      pause(0.00001)
      M(ct)=getframe(gcf);
      ct=ct+1;
end

%animation

movie(M)
videofile=VideoWriter(\'Simple Harmonic Motion.avi\',\'Uncompressed AVI\')
open(videofile)
writeVideo(videofile,M)
close(videofile)

The function which does the integration and gives a set of angular displ. and angular velocity:

function [dtheta_dt]= ode_fun(t,theta,k,g,l,m)
theta1=theta(1);
theta2=theta(2);
dtheta1_dt=theta2;
dtheta2_dt = -(k/m)*theta2-(g/l)*sin(theta1);
dtheta_dt=[dtheta1_dt;dtheta2_dt];
end

Here is the solution of the ODE using Matlab(Linear Algebra):

CODE:

clear all
clc
l=1.0;
b=0.5;
m=1.0;
g=9.81;
A=[-b/m -g/l;1.0 0.0]
[evec,eval]=eig(A)
initial=[3 0]'; % for t=0 theta_dot(0)=3, theta =0
t=linspace(0,20,600); % t is for 0 to 20 sec in the interval 0.1
Y=[evec(1) evec(3);evec(2) evec(4)]  %matrix with colums as eigen vector
c=inv(Y)*initial  % finding c1 and c2 using initial value
%r= Y*c
for i=t

   theta_dot=c(1)*exp(eval(1)*t)*evec(1)+ c(2)*exp(eval(4)*t)*evec(3);
   theta=c(1)*exp(eval(1)*t)*evec(2)+ c(2)*exp(eval(4)*t)*evec(4);
end
figure(1)
plot(t,theta_dot,'red')
hold on;
plot(t,theta,'blue')
ylabel('Responses')
xlabel('Time')
legend('velocity vs time','displacement vs time')
title('Using Matlab')

figure(2)
ct=1;
for i=1:length(t)
       THETA = theta(i);
    %Coordinates for pendulum
    x1 = l*sin(THETA);
    y1 = -l*cos(THETA);

     figure(2)
     % Plot string
     plot([0 x1],[0 y1],'linewidth',5)
     
     % Plot Support
     rectangle('Position',[-0.5 0 1 0.5],'EdgeColor','k','LineWidth',3)
     
     %plotting of the ball
     ball=viscircles([x1,y1],0.2);
     grid on
     
     
      axis([-2 2 -2 1])
      title('Simple Harmonic Motion in Pendulum Using Matrix')
      pause(0.00001)
      M(ct)=getframe(gcf);
      ct=ct+1;
end
%animation
movie(M)
videofile=VideoWriter('Simple Harmonic Motion1.avi','Uncompressed AVI')
open(videofile)
writeVideo(videofile,M)
close(videofile)

RESULTS AND CONCLUSION:

NUMERICAL INTEGRATION

USING MATRIX FORM

The plot describes that as time increases the velocity and displacement of the SHM decreases since we have used the damping factor. The fun fact is that the solution obtained from calculus and linear algebra are exactly identical for this Differential Equation system.

The link to the simulations are provided below:

Numerical Integration with c = 0.05:

Linear Algebra solution using c = 0.5: