Question 10.8: Damped, nonlinear motion of a pendulum. A pendulum is modele......
Damped, nonlinear motion of a pendulum.
A pendulum is modeled by a mass that is attached to a weightless rigid rod. According to Newton’s second law, as the pendulum swings back and forth, the sum of the forces that are acting on the mass equals the mass times acceleration (see the free body diagram in the figure). Writing the equilibrium equation in the tangential direction gives:
\Sigma F_{t}=-c L \frac{d \theta}{d t}-m g \sin \theta=m L \frac{d^{2} \theta}{d t^{2}} (10.153)
where \theta is the angle of the pendulum (with respect to the vertical axis, as shown in the figure), c=0.16(\mathrm{~N} \cdot \mathrm{s}) / \mathrm{m} is the damping coefficient, m=0.5 \mathrm{~kg} is the mass, L=1.2 \mathrm{~m} is the length, and g=9.81 \mathrm{~m} / \mathrm{s}^{2} is the acceleration due to gravity.
Equation (10.153) can be rewritten as the following second-order differential equation:
\frac{d^{2} \theta}{d t^{2}}=-\frac{c}{m} \frac{d \theta}{d t}-\frac{g}{L} \sin \theta (10.154)
The pendulum is initially displaced such that \theta=90^{\circ}, and then at t=0 it is released from rest, \frac{d \theta}{d t}=0 (zero initial velocity). Determine the angle of the pendulum as a function of time, \theta(t), for the first 18 seconds after it is released.
To solve Eq. (10.154), a new dependent variable, w, is introduced, such that:
w=\frac{d \theta}{d t} \quad \text { and } \quad \frac{d w}{d t}=\frac{d^{2} \theta}{d t^{2}}
With these definitions Eq. (10.154) can be rewritten as the following system of two first-order ODEs:
\begin{gathered}\frac{d \theta}{d t}=w \text { with the initial condition } \theta(0)=\frac{\pi}{2} & (10.155) \\\frac{d w}{d t}=-\frac{c}{m} w-\frac{g}{L} \sin \theta \text { with the initial condition } w(0)=0 & (10.156)\end{gathered}
Script File
The system of ODEs is solved with a user-defined MATLAB function named [t, x, y]= SYS2ODESRK4 (ODE 1, ODE 2, a, b, h, x 1, y 1 ), which is listed below. The function solves a system of two ODEs using the fourth-order Runge-Kutta method. The function uses two additional userdefined functions that calculate the value of the right-hand side of each ODE.
Program 10-6: User-defined function. Solving a system of two first-order ODEs using fourth-order Runge-Kutta method.
function [t, x, y]=Sys2ODEsRK4(ODE1,ODE2,a,b,h,x1,y1)
% Sys2ODEsRK4 solves a system of two first-order initial value ODEs usinq
% fourth-order Runqe-Kutta method.
% The independent variable is t, and the dependent variables are x and y.
% Input variables:
% ODE1 Name for the function that calculates dx/dt.
% ODE2 Name for the function that calculates dy/dt.
% a The first value of t.
% b The last value of t.
% h The size of a increment.
% x1 The initial value of x.
% y1 The initial value of y.
% Output variables:
% t A vector with the t coordinate of the solution points.
% x A vector with the x coordinate of the solution points.
% y A vector with the y coordinate of the solution points.
t(1) = a; x(1) = x1; y(1) = y1;
n = (b - a)/h;
for i = 1:n
t(i+1) = t(i) + h;
tm = t(i) + h/2;
Kx1 = ODEl(t(i) ,x(i) ,y(i));
Ky1 = ODE2(t(i) ,x(i),y(i));
Kx2 = ODE1(tm,x(i)+ Kx1*h/2,y(i)+ Ky1*h/2);
Ky2 = ODE2(tm,x(i)+ Kx1*h/2,y(i)+ Ky1*h/2);
Kx3 = ODE1(tm,x(i)+ Kx2*h/2,y(i)+ Ky2*h/2);
Ky3 = ODE2(tm,x(i)+ Kx2*h/2,y(i)+ Ky2*h/2);
Kx4 = ODE1(t(i + 1),x(i)+ Kx3*h,y(i)+ Ky3*h);
Ky4 = ODE2(t(i + 1),x(i)+ Kx3*h,y(i)+ Ky3*h);
x(i+1) = x(i) + (Kx1 + 2*Kx2 + 2*Kx3 + Kx4)*h/6;
y(i+1) = y(i) + (Ky1 + 2*Ky2 + 2*Ky3 + Ky4)*h/6;
end
To solve the system of ODEs in Eqs. (10.155) and (10.156), two user-defined functions are written. One, named PendulumDthethaDt, calculates \frac{d \theta}{d t}, and the other, named PendulumDthethaDt, calculates \frac{d w}{d t} :
function dxdt = PendulumDthethaDt(t,x,y)
dxdt = y;
function dydt = PendulumDwDt(t,x,y)
c = 0.16; m = 0.5; g = 9.81; L = 1.2;
dydt = -(c/m)*y - (g/L)*sin(x);
The user-defined functions are used in the following MATLAB program (script file) for solving the pendulum problem. The program also displays a plot ( \theta versus time) of the solution.
plot(t,x)
xlabel('Time (s) ')
ylabel ( 'Angle thetha (rad) ' )
When the script file is executed, the following figure is displayed:
