Week 3 - Solving second order ODEs
INTRODUCTION:
The presented project report attempts to provide a comprehensive study and develop a script for solving second order differential equations for the transient motion using MATLAB is software featuring A high - level programmin g language, primarily intended for numerical computations. octave helps in solving linear and nonlinear problems numerically, and for performing other numerical experiments to describe the transient behaviour of a system.
A simple example is a pendulum. the way the pendulum moves depends on the newton's second law. when this law is written down, we get a second order ordinary differential equation that describes of the "ball"w.r.t time.
THE SIMPLE PENDULUM:
A simple pendulum consists of a particle of mass attached to one end of a light inextensible string of length where the other end of the string is attached to a fixed point.The particle when at rest hangs vertically below the fixed point. the particle and string when displaced from the equilibrium position oscillate in a circular arc in the saame vertical plane. This physical description suggests the particle moves in 2D and therefore the motion will not behave like simple harmonic motion.
We consider a simple penduluym that consists of a mass suspended by a string or rod of length L (and negligible mass). The angle it makes with the vertical varies with time as a sine or cosine.
looking at the forces on the pendulum bob, we see that the restoring force is proportional to sin(theta),where as the restoring force for a spring is proportional to the displacement (which is theta in this case).
EXPLANATION TO THE CODE:
1) Initially a function (ode_pendulum) is created to the call the ordinally differential equation of the simple pendulum. the code for function.
function [dtheta_dt] = ode_pendulum(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
2) In the mainb program, i nitially the inputs such as damping coefficient ,length of the string, mass of the bob and acceleration due to gravity are decleared.
3)Then initial conditions of angular displacement and angular velocity are also given.
4)Time period is declared using the 'linspace' function to get multiple time points.
5) The ode_pendulum function is then called using ode45 a matlab solver for ode functions.
6) Now under the for loop we define the cooridinates (x0,y0)and(x1,y1) using the equations.
7)The plot command is used to plot the links on the graph. And to a animate the plot ,getframe command is used (getframe) returns a movieframe. and the counter variable is incrememented.
8) Both vfor loops are then ended and the syntax for movie is written.
MATLAB CODE FOR SIMPLE PENDULUM:
function h = pendulumMotion
% This function was created by saveas(...,'mfig'), or print -dmfile.
% It loads an HG object saved in binary format in the FIG-file of the
% same name. NOTE: if you want to see the old-style MATLAB code
% representation of a saved object, previously created by print -dmfile,
% you can obtain this by using saveas(...,'mmat'). But be advised that the
% MATLAB file/MAT-file format does not preserve some important information due
% to limitations in that format, including ApplicationData stored using
% setappdata. Also, references to handles stored in UserData or Application-
% Data will no longer be valid if saved in the MATLAB file/MAT-file format.
% ApplicationData and stored handles (excluding integer figure handles)
% are both correctly saved in FIG-files.
%
%load the saved object
[path, name] = fileparts(mfilename('fullpath'));
figname = fullfile(path, [name '.fig']);
if (exist(figname,'file')), h1=open(figname); else h1=open([name '.fig']);
end
if nargout > 0, h = h1;
end
SYTAX FOR THE MOVIE:
movie(M) %plays the movie in matrix M once.
videofile=videowritten ('videoname.avi',filetype) %assigns file name and file format
open (videofile,M)%writes the video with plot frames in the arrays M
close(videofile)%closes the videofile
MATLAB VIDEO CODING:
% Caleb Williams
% Animates the motion of a simple pendulum as a function of time.
clear all; close all; clc;
% All units in the SI system
m=1; % pendulum bob mass
l=2; % pendulum rod length
g=9.81; % acceleration due to gravity
t=10; % simulation time
c=0.25; % damping coefficient
tStep=0.01; % simulation time step
iterations=t/tStep; % total number of iterations
% pre-allocation of variables for speed
theta=zeros(iterations,1);
omega=zeros(iterations,1);
alpha=zeros(iterations,1);
T=zeros(iterations,1);
theta(1,:)=pi/2; % initial angular position
omega(1,:)=0; % initial angular velocity
alpha(1,:)=0; % initial angular acceleration
% Set up video
v=VideoWriter('pendulum.avi');
v.FrameRate=60;
open(v);
for n=1:iterations
% Iteratively solve equations of motion using Euler's Method
theta(n+1,:)=theta(n,:)+omega(n,:)*tStep; % new angular position
omega(n+1,:)=omega(n,:)+alpha(n,:)*tStep; % new angular velocity
alpha(n+1,:)=(-g*sin(theta(n+1,:)))/l-c*omega(n+1,:); % new angular acceleration
% Plot everything for the video
hold on;
fill([-l-0.2*l l+0.2*l l+0.2*l -l-0.2*l],[-l-0.2*l -l-0.2*l l+0.2*l l+0.2*l],'w'); % clears background
plot([0 l*sin(theta(n,:))],[0 -l*cos(theta(n,:))],'b','LineWidth',3); % plots rod
plot(l*sin(theta(n,:)),-l*cos(theta(n,:)),'Marker','o','MarkerSize',20,'MarkerFaceColor','r','MarkerEdgeColor','r'); % plots bob
%plot(0,0,'Marker','0','MarkerSize',5,'MarkerFaceColor','k','MarkerEdgeColor','k'); % plots pin
xlim([-l-0.2*l l+0.2*l]);
ylim([-l-0.2*l l+0.2*l]);
title('Simple Pendulum by Caleb Williams');
frame=getframe(gcf);
writeVideo(v,frame);
end
close(v);