Challenge

Week 4.1 - Solving second order ODEs

GGaneshram Devakumar9 Aug 2024

Objective - Writing a program in MATLAB that will simulate the pendulum motion for given span of time and given angular acceleration and velocity, with these results of angular velocity of pendulum over the time, we need to animate the pendulum motion.

Simple Pendulum -

Equation of motion is:

In the above equation,

b = damping coefficient.
g = gravity in m/s^2.
L = length of the pendulum in the meter.
m = mass of the ball in kg.

ODE Function -

Input -

b = 0.05; %% Damping Coefficient
g = 9.81; %% Gravitational Force in m/s^2
l = 1; %% Length of pendulum in m
m = 0.1; %% Mass in Kg
Start angular displacement in rad.
Start angular velocity in rad/sec^2
Time of oscillation - 0 to 20 secs

MATLAB Code -

%% Solving a Second Order Differential Equation For Simple Pendulum

clear all

close all

clc

% Input

b = 0.05; %% Damping Coefficient

g = 9.81; %% Gravitational Force in m/s^2

l = 1; %% Lenght of pendulam in m

m = 0.1; %% Mass in Kg

% Intial Condition

theta_0 = [0;3]; % [Start angular Displacement;Start angular Velocity]

% Time points

t_span = linspace(0,20,100); %% Time of oscillation between 0 to 20 sec with 100 milli sec

%solve ODE

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

plot(t,results(:,1),'LineStyle','--')

hold on

plot(t,results(:,2),'LineStyle',':')

xlabel('Time in sec')

ylabel('Plot')

title(['Second Order Differential Equation for simple Pendulum'])

legend('Angular Postion in Rad', 'Angular Velcity in rad/s^2')

% Creation of animation video of simple pendulam

ct=1 % Intializing the frame count as 1

theta1 = results(:,1);

for i = 1:length(theta1);

x0=0;

y0=0;

theta2 = theta1(i)

x1=1*sin(theta2);

y1=-1*cos(theta2);

% Poltting

plot([x0 x1], [y0,y1],'linewidth',4)

% string of pendulum

line([x0 x1],[y0 y1],'linewidth',3,'color','g')

hold on

% bob of pendulum

plot(x1,y1,'marker','o','markersize',20,'markerfacecolor','m')

axis([-1 1 -1.5 0.5])

pause(0.03) % pausing the results to 0.3 sec

hold off

M(ct)=getframe(gcf) %% calling out the figure to create the movie

ct=ct+1 %% to display the frame (ct) by frame in one graph and incrementing for all the position of bob

end

movie(M)

videofile=VideoWriter('Challenge_4_1.m', 'Uncompressed AVI');

open(videofile);

writeVideo(videofile,M);

close(videofile);Output #1

Output #1

Output #2

Animation Video -

Video Link - https://youtu.be/5o_ktdF69C0

Steps Followed in this Program -

Initialization of start condition of a simple pendulum and pendulum input parameters(b,g,l,m,Time)
To solve the ODE in which MATLAB inbuilt function (ode45) is used with theta, t, and ode_func - another function is created as a separate own function.
Results are plotted w.r.t time using plot command with proper legends, label, Plot Titles
By using the results of angular position (100 position of pendulum in a time span of 3 secs), for loop used to access each value of pendulum position plotted
Using ‘Line’ command string position plotted and Bob plotted with plot command
Position of pendulum (theta1 and theta2) for each angular position is created with for loop and increasing the ‘ct’ one by one for 100 times in a span of 3 sec
Then move command is used to create an animated video of the oscillating pendulum.

Error Faced in this Program -

#1 Pendulum oscillates without fixed position and oscillates in horizontal direction as shown below -

Direction of sin theta and cos theta (x1=1*sin(theta2);y1=-1*cos(theta2)) changed as per pendulum concept and window size 'axis([-1 1 -1.5 0.5])' added in a command to visualize the complete pendulum with fixed position

#2 Pendulum oscillates vibration looks small in a animation video -

The intimal start angular position (what we given) and final angular position(what we got from solving ODE function) of pendulum is in radians not in Degree. When assigning postion of pendulum wrongly given as theta " x1=1*sind(theta2); y1=-1*cosd(theta2);" below changes done in program

below changes done in program to solve this error x1=1*sin(theta2);y1=-1*cos(theta2);

Week 4.1 - Solving second order ODEs

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.