Solving second order ODEs by MATLAB

Objective

write a program that solves the following ODE. This ODE represents the equation of motion of a simple pendulum with damping.

Theory

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.

CODE

function

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

Main program

close all
clear all
clc

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

%initial condition
theta_0=[0 3];

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

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

plot(t,results(:,1))
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')

%animation of transient motion
figure(2)
ct=1;
x0=0;
y0=0;
outcome=results(:,1);

for i=1:length(outcome)
   
    pos=outcome(i); 
    x1=l*sin(pos);
    y1=-l*cos(pos);
    
    plot([-1.5 1.5],[0 0])
    hold on
    plot([x0 x1],[y0 y1],x1,y1,'g.','markersize',20,'linewidth',1,'color','b')
    axis([-1.5 1.5 -1.5 0.5])
    hold off
    M(ct)=getframe(gcf);
    ct=ct+1;
    pause(0.03)
  
  
end

%making the video file

movie(M)
videofile=VideoWriter('simple_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

CODE EXPLANATION

close all
clear all
clc

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

%initial condition
theta_0=[0;5];

%time points
t_span=linspace(0,10,500);

%solve ode
[t,results]=ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0); %function defination,ode45 is inbuilt passing t and theta main

plot(t,results(:,1)) %argument ode_func is calling and passing arguments.
hold on
plot(t,results(:,2))
ylabel('plot')
xlabel('time')

%animation of transient motion %figure 2 shownfigure(2)
ct=1; % counter set as 1
x0=0; %start of length set as origin
y0=0;
outcome=results(:,1); %storing result as outcome

for i=1:length(outcome) %for loop from 1 to 20 sec at end of length

pos=outcome(i);
x1=l*sin(pos); %end point of length
y1=-l*cos(pos);

plot([-1.5 1.5],[0 0]) %plotiing
hold on
plot([x0 x1],[y0 y1],x1,y1,'g.','markersize',20,'linewidth',1,'color','b')
axis([-1.5 1.5 -1.5 0.5])
hold off
M(ct)=getframe(gcf);
ct=ct+1;
pause(0.03)


end

%making the video file

movie(M)
videofile=VideoWriter('simple_pendulum.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)

OUTPUT

Animation youtube link

https://youtu.be/1X5RrWJb8PE