Solving 2nd Order Ordinary Differentiation Equation using MATLAB

A program to simulate the transient behaviour of a simple pendulum and to create an animation of it's motion

ode45 is a versatile ODE solver and is the first solver you should try for most problems. However, if the problem is stiff or requires high accuracy, then there are other ODE solvers that might be better suited to the problem.

A function is a group of statements that together perform a task. Functions operate on variables within their own workspace, which is also called the local workspace, separate from the workspace you access at the MATLAB command prompt which is called the base workspace.

MAIN MOTIVE :

-Use ODE 45 to solve the equation

-Insert Function as data for the equation

-Find the relationship between Time and Displacement

-Simulate the motion of pendulum as per data calculated by the MATLAB

PROCEDURE :

-Basic commands to clear,close,clear command window of previous data

-Input values of Damper,Force of gravity,Length of pendulum,Mass of Pendulum

-Input Initial Condition,

i.e,

Any variable = [Displacement(can be known as theta1);Velocity(can be known as Velocity)]
theta0=[0;15];

Limits for integration,
time_limits=linspace(0,100,500);

-CALCULATION OF ORDINARY DIFFERENTIATION EQUATION(Application of ODE 45, and calling Function)
[t(time) (independant variable),results=(theta1+theta2)(dependant variable)]
[t,results]=ode45(@(t,theta)ode_function(t,theta,d,g,l,m),time_limits,theta0)

%PLOTTING RESULTS IN GRAPH
subplot(2,1,1) - Divides the graph sheet to add another figure in same window to relate the data easy
figure(1)
plot(t,results(:,1)) - Plot time data
hold on
plot(t,results(:,2))-Plot theta 1 column from the result matrix which is combination of theta 1 & theta 2
xlabel("Time")
ylabel("Displacement")
plot(t,results(:,1),'o','markerfacecolor','black','MarkerSize',2) - Adds points along the plots
plot(t,results(:,2),'o','markerfacecolor','black','MarkerSize',2) - Adds points along the plots
hold off

ct=1;

ANIMATING THE MULTIPLE GRAPH PLOT IN SINGLE SCREEN
%Simple pendulum
for i=1:1:length(time_limits);
A=results(i,1);
x=0;
y=0;
x1=l*sind(A);
y1=-l*cosd(A);

%Plotting & Simulating Pendulum in graph

subplot(2,1,2)
plot([x,x1],[y,y1],'linewidth',5,'color','r')
hold on
plot(x1,y1,'o','markerfacecolor','black','MarkerSize',25)
plot(x,y,':bs','markerfacecolor','black','MarkerSize',50)
hold off
axis([-1 2 -1.5 .1])%Graph data visibility for animation or simply axis limit
pause(0.05)%pausing between each graph in which speed at which the graph animation is to be met
Z(ct)=getframe(gcf);% to get graph data at current
ct = ct+1;
end

%MOVIE COMMAND
movie(Z)
Videofile=VideoWriter('ode.avi','uncompressed AVI');
open(Videofile)
writeVideo(Videofile,Z)
close(Videofile)

MAIN PROGRAM

clear all
close all
clc


%INPUTS
 %Damper,
 d=0.05;
 %Force of gravity
 g=9.81;
 %Length of pendulum,
 l=1;
 %Mass of Pendulum,
 m=1;
 %Initial Condition,
 %Any variable = [Displacement(can be known as theta1);Velocity(can be known as Velocity)]
 theta0=[0;15];
 %Limits for integration,
 time_limits=linspace(0,100,500);
 
%CALCULATION OF ORDINARY DIFFERENTIATION EQUATION:
%[t(time) (independant variable),results=(theta1+theta2)(dependant variable)]
[t,results]=ode45(@(t,theta)ode_function(t,theta,d,g,l,m),time_limits,theta0)


%PLOTTING RESULTS IN GRAPH
subplot(2,1,1)
figure(1)
plot(t,results(:,1))
hold on
plot(t,results(:,2))
xlabel("Time")
ylabel("Displacement")
plot(t,results(:,1),'o','markerfacecolor','black','MarkerSize',2)
plot(t,results(:,2),'o','markerfacecolor','black','MarkerSize',2)
hold off

ct=1;


%Simple pendulum
for i=1:1:length(time_limits);
    A=results(i,1);
    x=0;
    y=0;
    x1=l*sind(A);
    y1=-l*cosd(A);
    

    
%Plotting & Simulating Pendulum in graph

subplot(2,1,2)
plot([x,x1],[y,y1],'linewidth',5,'color','r')
hold on
plot(x1,y1,'o','markerfacecolor','black','MarkerSize',25)
plot(x,y,':bs','markerfacecolor','black','MarkerSize',50)
hold off
axis([-1 2 -1.5 .1])%Graph data visibility for animation or simply axis limit
pause(0.05)%pausing between each graph in which speed at which the graph animation is to be met
Z(ct)=getframe(gcf);% to get graph data at current
ct = ct+1;
end


%MOVIE COMMAND
movie(Z)
Videofile=VideoWriter('ode.avi','uncompressed AVI');
open(Videofile)
writeVideo(Videofile,Z)
close(Videofile)


    

FUNCTIONAL PROGRAM :

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

ERROR OCCURED & SOLVED :

RESULTS :

Simulation Link : Pendulum simulation https://youtu.be/YYhyQkqf49Y