Animaciones con Matlab/GNU-Octave

Esto puede resultar de ayuda cuando se quiere calcular la dinámica de un sistema mecánico, es decir su posición en función del tiempo, su velocidad o su aceleración para el cálculo de las fuerzas de trepidación. Además de el cálculo, uno puede "ver" estas magnitudes, en un animación.

Este es el ejemplo clásico del cálculo de un mecanismo de cuatro barras. Se llama de cuatro barras, pero sólo tres de ellas se mueven

Animación en Matlab o GNU/Octave

Simplemente es ir calculando para cada instante la posición de nuestro "sistema", las guardamos en matrices o vectores, y después graficamos para cada instante. Luego si se desea se puede guardar en un gif y se ve como el pistón de arriba ( para ver más Aquí para hacer un gif con matlab

Conocer la posición de un sistema Pistón Biela Manivela

Se plantea la posición del sistema en un instante cualquiera para establecer las ecuaciones que permitan calcularla.

La geometría conocida del sistema es:

• x posición del pistón respecto del eje de giro de la manivela


• l_b Longitud de la biela


• l_m Longitud de la manivela


• t₁ Ángulo entre la línea de recorrido del pistón y la manivela


• t₂ Ángulo entre la línea de recorrido del pistón y la biela

Las ecuaciones se determinan sumando en cada componente los desplazamientos a partir de un punto y llegando al mismo punto. En este caso, partiendo del punto de rotación de la manivela, los desplazamientos serían:

Es sistema de ecuaciones queda:




Conocida la geometría, el largo de la biela y la manivela, y dado el ángulo de la línea de recorrido del pistón y la manivela, el sistema de ecuaciones que queda, tiene como incógnita la posición del pistón y el ángulo.


Se puede resolver despejando en la primera ecuación el valor de t₂ y luego reemplazar en la segunda ecuación para encontrar x. Pero para generalizar el método, es conveniente resolver de otra forma el sistema para casos en que hay dos incógnitas, ambas en las dos ecuaciones y como argumento de una función trigonométrica, como en el caso del mecanismo de cuatro barras.

Resolución del sistema de ecuaciones

Para ello se define una función F de dimensión 2, que tiene componentes los miembros del sistema de ecuaciones distintos de cero:


La función se iguala a cero por el segundo miembro del sistema de ecuaciones.
Para encontrar las incógnitas se utiliza la aproximación de la función por los dos primeros términos de la serie de Taylor.



Lo que desarrollado es:



Se impone que el valor de la aproximación de la serie de Taylor sea cero, coincidente con el de la función, entonces queda:


Escrito de otra manera:


Por lo que basta con invertir la matriz para conocer el valor de las incógnitas de la aproximación de la función por la serie de Taylor. Dichos valores son una aproximación a los valores que reemplazados en el sistema lo satisfacen. Por lo que se programa una iteración, de manera de partir con ciertos valores iniciales, evaluar la función, calcular las derivadas parciales, invertir la matriz, obtener valores nuevos, sumarlos a los antiguos, evaluar la función para que sea menor a un error pequeño predefinido y si es así estos valores son considerados los que satisfacen el sistema. Si no es así, se vuelve a iterar. Lo que aparece resumido en el diagrama de flujo.

La función de cálculo

Es diagrama de flujo anterior implementado en código GNU-Octave/Matlab se muestra a continuación. Nótese que como valores de entrada están el ángulo de la manivela, para el cual el sistema tiene una única posición y además se incluyen los valores de x y t₂ iniciales, los que deben ser tentativos.
La función entrega como parametro de salidas las incógnitas x y t₂ y ademas de esos dos vectores, xx e yy, en los que se guardan la posición de los puntos para cada instante en los ejes x e y.

function [xx,yy,x,t2] = blogpiston(t1,x,t2) % Tolerancia del error e=0.001; % Dimensiones del sistema lm=35; lb=56; % definicion de F f1=lm*sin(t1)-lb*sin(t2); f2=lm*cos(t1)+lb*cos(t2)-x; while abs(f1)>e | abs(f2)>e df1x=0; df1t2=-lb*cos(t2); df2x=-1; df2t2=-lb*sin(t2); A=[df1x,df1t2;df2x,df2t2]; B=[-f1;-f2]; % Resolucion del sistema X=A\B; % Por si acaso no converge if X==[0;0] break end %valores nuevos de x y t2 x=x+X(1); t2=t2+X(2); %definicion de F con los nuevos x y t2 f1=lm*sin(t1)-lb*sin(t2); f2=lm*cos(t1)+lb*cos(t2)-x; end % Puntos de la maya para el t1 dado xx=[0,lm*sin(t1),0]; yy=[0,lm*cos(t1),x];

La animación

Para animar el sistema debemos tener calculada la posición del sistema dados varios instantes, o mejor dicho, para varios ángulos de la manivela. Así, si se supone una velocidad angular constante en la manivela, puede calcularse para ángulos t₁ conocidos, a un paso conocido,y luego graficar cada resultado del sistema para cada valor de angulo t₁.

Las instrucciones anteriores se indican en el siguiente código:

% animacion clear all; % paso paso=6*pi/180; % numero de puntos nc=100; % datos t1=pi/4; % estimaciones x=70; t2=30*pi/180; % definicion de los vecotres para guardar los puntos X=[0,0,0]; Y=X; XP=[0,0,0,0,0]; YP=XP; % calculo para cada posicion de t1 for i=1:nc [xx,yy,x,t2] = blogpiston(t1,x,t2); X(i,:)=xx; Y(i,:)=yy; t1=t1+paso; %cuadrado para simular piston (Es para que se vea bonito ) [xx,yy]=cuad(xx(3),yy(3),5,0); XP(i,:)=xx; YP(i,:)=yy; % end % grafico/animacion for i=1:nc %Marco para la pantalla plot(-100,-100) hold on plot(100,100) %malla plot(X(i,:),Y(i,:),'-o') %cuadrado piston plot(XP(i,:),YP(i,:)) %tiempode pausa pause(.1) hold off end

Lo anterior puede guardarse como imágenes y puede elaborarse un gif de un pistón agregando la sintaxis del artículo citado, y quedará una animación como la que se ve en la cabeza del artículo.

Dibujar el pistón

Como la "animación" descrita dibuja solo lineas y puntos, agregué una función que dibuja un cuadrado para representar el pistón, tomando como referencia el punto calculado para la posición de este

Notar que x e y son las coordenadas del pistón. El pistón es un cuadrado de lado e que recorre una línea horientada con la vertical de ángulo a

function [X,Y]=cuad(x,y,e,a); r=e*sqrt(2); b=pi/4-a; X=[x-r*sin(b),x+r*cos(b),x+r*sin(b),x-r*cos(b),x-r*sin(b)]; Y=[y-r*cos(b),y-r*sin(b),y+r*cos(b),y+r*sin(b),y-r*cos(b)];

Archivos *.m

Dada cualquier eventualidad los archivos *.m están en los siguientes enlaces:

• Función


• Programa


• Cuadrado para el pistón

Francisco Roco