Week 3 - Solving second order ODEs
Aim:
To Write a program in Matlab\Octave that will simulate the pendulum motion & write a program that solves the following ODE.
Simple Pendulum:
A pendulum is a body suspended from a fixed support so that it swings freely back and forth under the influence of gravity. When a pendulum is displaced sideways from its resting, equilibrium position, it is subject to a restoring force due to gravity that will accelerate it back toward the equilibrium position. When released, the restoring force acting on the pendulum's mass causes it to oscillate about the equilibrium position, swinging it back and forth. The mathematics of pendulums are in general quite complicated. Simplifying assumptions can be made, which in the case of a simple pendulum allow the equations of motion to be solved analytically for small-angle oscillations.
The differential equation which represents the motion of a simple pendulum is
`(d^2theta)/dt^2 +( b/m)(d(theta))/(dt) + g/Lsintheta = 0`
where ; b= damping coeficient
m= mass of bob
g= acceleration due to gravity
L= length of the rod connecting bob
Assuptions- :
Rod is massless
bob is a point mass
Main Program explanation-:
These equations are used in deriving pendulum motion.
- `(d(theta1))/dt =theta2`
2. `(d^2theta)/dt^2 = (d(theta2))/dt`
3. `(d(theta2))/dt= -(b/m).theta2- (g/L)sintheta1`
4. `d/dt[theta1 ,theta2] = [theta2 ,-(b/m)theta2- (g/L)sintheta1]`
4. is a solution vector
- In the main program, initially the inputs such as damping coefficient, length of the string , mass of the bob and acceleration due to gravity are declared.
- Then initial conditions of angular displacement and angular velocity are also given.
- Time period is declared using the 'linspace' function to get multiple time points.
- ODE means ordinary differential equation
- The ode_pendulum function is then called using ode45 a matlab solver for ode functions
-
ode45
Solve nonstiff differential equations — medium order method,
ode45is a versatile ODE solver and is the first solver you should try for most problems. However, if the problem is stiff or requires high accuracy, then there are other ODE solvers that might be better suited to the problemwe have to provide two values for theta_0 i,e [ 0; 5] for two ode means for 1st order ode and 2nd order ode respectively.
- Now under the for loop we define the coordinates (X0, Y0) and (X1, Y2) using the equations.
- The plot command is used to plot the links on the graph. And to animate the plot, getframe command is used (getframe returns a movie frame. The frame is a snapshot (pixmap) of the current axes or figure). And the counter variable is incremented.
- Both for loops are then ended and the syntax for movie is written.
MAIN PROGRAM
clear all
close all
clc
%inputs
b= 0.02;
g= 9.81;
l= 1;
m= 0.1;
%inital condition
%theta_0= [disp;velo]
theta_0= [0; 5]; % theta_0= [initial displace ment; initial velocity]
%time points
t_span= linspace(0,10,500); %time period for different time points
%solve ode, syntax used for ode45
[t, results]= ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span, theta_0);
plot(t, results(:,1)) %plot b/w time & displ
hold on
plot(t,results(:,2)) %plot b/w time & displ
ylabel('plot')
xlabel('time')
ct= 1;
for i= 1:length(results) %loop to obtain diffr orientation of pendulum
figure(2) %calculatng coordinates of pendulum
x1= l*sin(results(i, 1));
y1= l*cos(results(i, 1)); % sin not sind
plot([-1, 1], [0, 0], [0, -x1], [0,-y1],'o','markers',15,'markerfacecolor','r')
line([0,-x1],[0,-y1],'linewidth',4)
axis([-2 2 -2 2]);
title('simple pendulum by rajat')
pause(0.01)
xlabel('Iteration');
ylabel('Tension (N)');
title('Tension in Pendulum Rod as a Function of Time');
M(ct)= getframe(gcf) %storing frames into M
ct= ct+1;
end
movie(M); %TO CREATE AVI FORMAT ANIMATION
videofile= VideoWriter('solving_2_nd_diffrential_equation_pendulum_motion.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile, M)
close(videofile)
Function Code
Initially a function (ode_func) is created to the call the ordinary differential equation of the simple pendulum saved in the array called as dtheta_dt. The code for 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
Results:
OUTPUT:
ERROR FACED:
dot indexing error, its normal error occurs when the file doesnt run completely
Referrence:
https://in.mathworks.com/help/matlab/ref/ode45.html
https://en.wikipedia.org/wiki/Pendulum
YOUTUBE'S LINK:
https://youtu.be/-kwFRPQItUA