Solving 2nd Order Ordinary Differential Equations using Matlab

Introduction

In Engineering, ODE is used to describe the transient behavior of a system. A simple example is a 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. Similarly, there are hundreds of ODEs that describe key phenomena and if you have the ability to solve these equations then you will be able to simulate any of these systems.

The formula for 2nd Order ODE

`(d^2θ)/(dt)^2+(b/m)*(dθ)/(dt)+(g/L)*sin(θ)=0`

In the above equation,

`g =`gravity in `m*s^-2`,

`L =`length of the pendulum in `m`,

m`=` mass of the ball in `kg`,

`b=`damping coefficient,

`θ=`Angle of the pendulum with horizontal.

Code

1. First input parameters are assigned and then the initial conditions.
2. A time matrix is also assigned from `0` to `20s` with `1000` points between them.
3. The time and results(angular position and velocity) are calculated using the ode45 function which takes input parameters such as the output of another function which calculates the differentiation of variable `θ` w.r.t time and time matrix and initial conditions.

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

4. Angular position is plotted alongside with angular velocity.
5. Angular position vs angular velocity is plotted.
6. Using for loop and angular position the motion of the pendulum is animated.
7. The animation is saved then.

% A program to solve 2nd order ODE 
close all
clear all
clc


% Input parameters
b = 0.05
g = 9.81
l = 1
m = 1


% initial conditions
theta_o = [0;3]
%size(theta_o)


% time points
time_span = linspace(0,20,1000);

% solve ODE
[t, Results] = ode45(@(t,theta) Pendulum_Damped_Code1(t,theta,b,m,g,l),time_span,theta_o);

subplot(5,4,[13,14,17,18]);
plot(t,Results(:,1),'color','k');
hold on;
grid on;
plot(t,Results(:,2),'color','r');
xlabel('Time (sec)');
ylabel('Position and Angular Velocity');
hold off;
subplot(5,4,[15,16,19,20]);
plot(Results(:,1),Results(:,2));
xlabel('Position-Radians');
ylabel('Angular Velocity');
grid on;


% To animate the pendulum
n = 1;

for i=1:length(Results(:,1))
  x0=0;
  y0=0;
  x1=l*sin(Results(i,1));
  y1=-l*cos(Results(i,1));
  subplot(5,4,[1,2,3,4,5,6,7,8,9,10,11,12]);
  plot([-l l],[0 0],'linewidth',4,'color','k');
  hold on;
  line([x0 x1],[y0 y1],'linewidth',2,'color','b');
  hold on;
  plot(x1,y1,'-o','markers',20,'markerfacecolor','r');
  
  grid on;
  axis([-1.5 1.5 -1.5 0.5]);
  pause(0.0003);
  hold off;
  q(n)=getframe(gcf);
  n=n+1;
  
endfor
movie(q)
videofile=VideoWriter('Pendulum_Damped.avi','uncompressed avi');
open(videofile)
writeVideo(videofile,q)
close(videofile)

Errors

The input parameters to calculate differentiations of θ were not the same and the plots were wrong.

Output

Input parameters and initial conditions are displayed:

Plots:

Animation with other plots:

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