octave :: single pendulum

octave_mechanics_singlependulum

\[T=\frac{1}{2}ml\dot{\theta}^2\] \[U=mgl(1-\cos\theta)\] \[L=T-U=\frac{1}{2}ml\dot{\theta}^2-mgl(1-\cos\theta)\] \[\frac{d}{dt}\frac{dL}{d\dot\theta}-\frac{dL}{d\theta}=0\] \[\ddot\theta =-\frac{g}{l}\sin\theta\] \[ \frac{d}{dt} \left[ \begin{array}{cc} \theta(t)\\ \dot\theta(t) \end{array} \right] = \left[ \begin{array}{cc} \dot\theta(t)\\ -\frac{g}{l}\sin\theta(t) \end{array} \right] , \left[ \begin{array}{cc} \theta(0)\\ \dot\theta(0) \end{array} \right] = \left[ \begin{array}{cc} \theta_0\\ 0 \end{array} \right] \]

Function: octave_mechanics_singlependulum ()

function out = octave_mechanics_singlependulum(m,l,g,x0,x1)
dt =  0.01;
%m = 1;
%l = 0.5;
%g = 9.8;
%x_k = [30/180*pi ;
%          0        ];
x_k = [x0 ;
       x1];
x = [];
tim =[];

for t = 0 : dt : 10
    k1 = [ x_k(2);
          -g/l * sin(x_k(1)) ] * dt;
    x_tmp = x_k + k1/2;
    k2 = [ x_tmp(2);
          -g/l * sin(x_tmp(1)) ] * dt;
    x_tmp = x_k + k2/2;
    k3 = [ x_tmp(2);
          -g/l * sin(x_tmp(1)) ] * dt;
    x_tmp = x_k + k3;
    k4 = [ x_tmp(2);
          -g/l * sin(x_tmp(1)) ] * dt;
    x_k1 = x_k + (k1 + 2 * k2 + 2 * k3 + k4 ) /6 ;
    x = [x x_k];
    tim = [tim t];
    x_k = x_k1;
endfor

plot(tim, x(1, :),tim, x(2, :));
out=0;
endfunction

octave_mechanics_singlependulum.m

octave_mechanics_singlependulum