Question 7.14: RK4 METHOD FOR SYSTEMS The pendulum system in Figure 7.5 con...
RK4 METHOD FOR SYSTEMS
The pendulum system in Figure 7.5 consists of a uniform thin rod of length l and a concentrated mass m at its tip. The friction at the pivot causes the system to be damped. When the angular displacement θ is not very small, the system is described by a nonlinear model in the form
\frac{3}{4}ml^{2} \overset{\cdot \cdot }{\theta }+0.18\overset{\cdot }{\theta } +\frac{1}{2}mgl\sin \theta =0Assume, in consistent physical units, that ml² = 1.28, \frac{g}{l}=7.45.
1. Transform the model into a system of first-order initial-value problems.
2. Write a MATLAB script that utilizes the user-defined function RK4System to solve the system in (1). Two sets of initial conditions are to be considered: (1) θ(0) = 15°, \overset{\cdot }{\theta } (0)=0 and (2) θ(0) = 30°, \overset{\cdot }{\theta } (0)=0 . The file must return the plots of the two angular displacements corresponding to the two sets of initial conditions versus 0 ≤ t ≤ 5 in the same graph. Angle measures must be converted to radians. Use at least 100 points for plotting purposes.
1. First note that
mgl=(ml²)(\frac{g}{l} )=(1.28)(7.45)=9.5360Therefore, the equation of motion is rewritten as
\frac{3}{4}(1.28)\overset{\cdot \cdot }{\pmb{\theta } } +0.18 \overset{ \cdot }{\pmb{\theta } }+\frac{1}{2}(9.5360)\sin \theta =0\Rightarrow 0.96\overset{\cdot \cdot }{\pmb{\theta } }+0.18\overset{ \cdot }{\pmb{\theta } }+4.7680\sin \theta =0There are two state variables: u_1 = \theta , u_2 = \overset{\cdot }{\pmb{\theta } }. The state-variable equations are formed as
\overset{\cdot }{u_1} =u_2
\overset{\cdot }{u_2}=-0.1875u_2-4.9667 \sin u_1
The two sets of initial conditions are
In vector form,
\overset{\cdot }{\textbf{u}}=\textbf{f}(t,\textbf{u}), u=\left\{\begin{matrix} u_1 \\ u_2 \end{matrix} \right\} , \textbf{f}(t,\textbf{u})=\left\{\begin{matrix} u_2 \\ -0.1875u_2-4.9667\sin u_1 \end{matrix} \right\}, \textbf{u}_0=\left\{\begin{matrix} \frac{15}{180}\pi \\ 0 \end{matrix} \right\} and \left\{\begin{matrix} \frac{30}{180}\pi \\ 0 \end{matrix} \right\}2.
t = linspace(0,5); u0 = [15*pi/180;0];
f = @(t,u)([u(2);−0.1875*u(2)−4.9667*sin(u(1))]);
uRK4 = RK4System(f,t,u0);
theta1 = uRK4(1,:); plot(t,theta1)
hold on
u0 = [30*pi/180;0];
uRK4 = RK4System(f,t,u0);
theta2 = uRK4(1,:); plot(t,theta2) % Figure 7.6
