Week 3 - Solving second order ODEs

second order ode

AIM = to solve the simple harmonic motion of simple pendulum ode

theory:

First we will employ our understanding of the underlying physics of the simple pendulum system to derive the structure of the system model. We will term this process, "modeling from first principles." In this example we employ the following variables.

(m)       mass of the pendulum bar

(M)       mass of the pendulum end weight

(l)       length to end weight center of mass

(theta)   pendulum angle from vertical (down)

To begin, we first draw the free-body diagram where the forces acting on the pendulum are its weight and the reaction at the rotational joint. We also include a moment due to the friction in the joint (and the rotary potentiometer). The simplest approach to modeling assumes the mass of the bar is negligible and that the entire mass of the pendulum is concentrated at the center of the end weight.

The equation of motion of the pendulum can then be derived by summing the moments. We will choose to sum the moments about the attachment point since that point is the point being rotated about and since the reaction force does not impart a moment about that point.

(1)

Assuming that the mass of the pendulum is concentrated at its end mass, the mass moment of inertia is . A more accurate approach would be to consider the rod and end mass explicitly. In that case, the weight of the system could be considered to be located at the system's mass center . In that case, the mass moment of inertia is . Depending on the parameters of your particular pendulum, you can assess if this added fidelity is necessary.

For the experiment we will perform shortly, the simple pendulum employed consists of a rod of length and mass with an end mass of . Therefore, the difference between and is significant enough to include. The difference between and is also significant enough to include.

We will also initially assume a viscous model of friction, that is, where is a constant. Such a model is nice because it is linear. We will assess the appropriateness of this model later. Sometimes the frictional moment is not linearly proportional to the angular velocity. Sometimes, the stiction in the joint is significant enough that it must be modeled too.

Taking into account the above assumptions, our equation of motion becomes the following.

PROCEDURE : first give the variabels and then inputs , initial condition after that solve the ode by using slove ode command ,write the function command and then call the function in main code and for loop the each value then create the animation and movie of the code

CODE:

% second order ode to solve the pendulam motion
close all 
clear all
clc 
b = 0.03;
g = 9.81;
l = 2;
m = 0.5;

% initial condition 
theta_0=[0;3];

% time points
t_span=linspace(0,25,100);

% solve ODE
[t,result]= ode45(@(t,theta)ode_func(t,theta,b,g,l,m),t_span,theta_0);

plot(t,result(:,1))
hold on 
plot(t,result(:,2))


ct=1;
for i=1:length(result(:,1))
    
    x0=0;
    y0=0;
    
    x1=l*sin(result(i,1));
    y1=l*cos(result(i,1));
    
    figure(2)
    %plotting
     plot([1 -1],[0 0],'linewidth',3)
    hold on
    
    plot([-x0 -x1],[-y0 -y1],'linewidth',4)
    axis([-4 4 -4 4])
    plot(-x1,-y1,'.','markersize',35)
    M(ct)=getframe(gcf);
    ct=ct+1;
    hold off
end
movie(M)
videofile=VideoWriter('ode_function_for_loop.avi','uncompressed AVI');
open(videofile)
writeVideo(videofile,M)
close(videofile)













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

RESULT : video link

https://youtu.be/-y8ToUkbNpg