Week 3 - Solving second order ODEs

Objective: The Objective of this challenge is to write a  program that can solve the second-order ODE of a simple pendulum and make the animation for the moment of the pendulum in terms of angular displacement.

Introduction:

Many physical problems from fluid mechanics, heat transfer, rigid body dynamics, and elasticity are modeled by second-order PDE's. In some problems fourth-order PDE's do arise, however, as we split higher order ordinary differential equations into a system of first-order equations, it is also a common practice to split a 4th order PDE into two-second order PDE's along with the necessary boundary and initial conditions and solve them together.

Problem statement; the given condition is we need to write a code for a pendulum and extract the video out of it by using the inbuikt matlab comands.

given parameters are;

length of the pendulum = l=1m

mass of the ball = m = 1Kg

damping co-efficient = b = 0.05

gravity = g = 9.81 m^2/sec

time = t = [0-20] sec

angular dispcement = 0

angular velocity = 3 rad/sec^2  (at t = 0)

given equation is

we need to do some assumptions and hand calculation to get the 'theta1' and 'theta2'

The given equation is `(d^2theta)/dt^2+b/m*(d(theta))/dt+g/L*sintheta`

Let's say

`theta = theta_1`

`(d(theta))/(dt) = theta_2`

`(d^2theta)/dt = (d^2theta_1)/dt^2 = d/dt (d(theta))/dt = d/dt(theta_2)`

ODE's in MATLAB:

The Ordinary Differential Equation (ODE) solvers in MATLAB^® solve initial value problems with a variety of properties. The solvers can work on stiff or nonstiff problems, problems with a mass matrix, differential algebraic equations (DAEs), or fully implicit problems.

An ordinary differential equation (ODE) contains one or more derivatives of a dependent variable, y, with respect to a single independent variable, t, usually referred to as time. The notation used here for representing derivatives of y with respect to t is y′ for a first derivative, y′′ for a second derivative, and so on. The order of the ODE is equal to the highest-order derivative of y that appears in the equation.

For example, this is a second-order ODE:

In an initial value problem, the ODE is solved by starting from an initial state. Using the initial condition, y0, as well as a period of time over which the answer is to be obtained, (t0,tf), the solution is obtained iteratively. At each step, the solver applies a particular algorithm to the results of previous steps. At the first such step, the initial condition provides the necessary information that allows the integration to proceed. The final result is that the ODE solver returns a vector of time steps t=[t0,t1,t2,...,tf] as well as the corresponding solution at each step y=[y0,y1,y2,...,yf].

Reference: https://www.mathworks.com/help/matlab/math/choose-an-ode-solver.html

How to create an animation using MATLAB:

To create a video from the Matlab inbuilt commands, we need to follow the process shown below. Make sure that we don't make any typo's and regarding upper and lower cases of the letters.

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

Main Code:

% code to simulate the pendulum motion

clear all
close all
clc

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


% we need to represent initial conditions to the problem
theta_0 = [0 3]

% given time interval or time span
time_given = linspace(0,20,400)

% write the ODE45 solver now

[t, outputs] = ode45(@(t, theta) second_order_ode_pendulum(t, theta, g,l,m,b), time_given,theta_0)

figure(1)
plot(t, outputs(:,1))
hold on
plot(t, outputs(:,2))
xlabel ('Time')
ylabel ('Angular displacement vs Angular velocity')
legend ('Angular displacement','Angular velocity')
grid on

figure(2)
plot(outputs(:,1), outputs(:,2))
xlabel ('position')
ylabel ('Angular velocity')
legend('Angular velocity')

theta2 = outputs(:,1)
ct = 1
 for i = 1:length(theta2)
     theta3 = theta2(i)
     % for fixed bar and movung bar of a pendulum
     x0 = 0
     y0 = 0
     
     x1 = 1*sin(theta3)
     y1 = -1*cos(theta3)
     % giving axis, color for the pendulum bob, etc.,
     figure(3)
     plot([-1.1 1.1],[0,0],'linewidth',3)
     hold on
     plot([x0 x1],[y0 y1], 'g','linewidth',3)
     hold on
     plot(x1,y1,'o','markers',20,'markerfacecolor','m')
     xlabel('x-axis')
     ylabel('y-axis')
     hold off
     grid on
     axis ([-1.5 1.5 -1.5 0.5])
     % time frame to get the the video
     pause(0.003)
     M(ct) = getframe(gcf)
     ct = ct+1
 end
 % code to run and save the video 
 movie(M)
 videofile = VideoWriter('Simple Pendulum.avi','uncompressed AVI')
 open(videofile)
 writeVideo(videofile,M)
 close(videofile)
 

Function Code:

% splitting the second order derivative into two prts
function [dtheta_dt] = second_order_ode_pendulum (t,theta,g,l,m,b)

theta1 = theta(1)

theta2 = theta(2)

% need to look the equation in two two parts
dtheta1_dt = theta2

dtheta2_dt = (-(b/m)*theta2)-((g/l)*sin(theta1))

dtheta_dt = [dtheta1_dt;dtheta2_dt]
end

Results:

Time vs Angular displacement and Angular velocity  :

• The plot explains how the angular velocity and angular displacement is varying with respect to time

• We can observe here how angular velocity and displacement are changing with respect to time.

Position vs Angular velocity:

• The plot explains to us the position of the bob or ball on the scale of [-1 1] while it's angular velocity is going to decrease gradually.
• The bob in the simple pendulum starts at a point at a velocity of 3 m/s and as the solution starts, the bob slowly lost its momentum and made the impression of the path (which is shown in above plot).
• After the given time of 20 seconds, the velocity reaches to nearly 1.3 m/s.

Animation: https://youtu.be/BZJxRoOlSdU

Error:

• The error raises because we are not using the parenthesis symbol (@) in front of the line where the second-order function code gets solved.