Oscilaciones


Inicio
 Oscilaciones
 Forzadas
 Tipos de fuerza

Respuesta de un oscilador
a una fuerza impulsiva (II)
En esta página, continuamos con el estudio del comportamiento de un oscilador bajo la
acción de un fuerza impulsiva (una fuerza que se aplica durante un determinado tiempo)

El ejemplo más significativo es el pulso media onda

Pulso triangular

Sea la fuerza

F(t)={F0at0≤t≤a0t>aF(t)={F0at 0≤t≤a0 t>a

Su transformada de Laplace es

f(t)=F(t)m=(u(t)−u(t−a))F0mta=f0a(t−u(t−a)(t−a)−u(t−a)a)g(s)=f0a(1s2−exp(−as)s2−aexp
(−as)s)f(t)=F(t)m=(u(t)−u(t−a))F0mta=f0a(t−u(t−a)(t−a)−u(t−a)a)g(s)=f0a(1s2−
exp(−as)s2−aexp(−as)s)

Para el caso de que las condiciones iniciales sean: posición inicial x0=0 y velocidad
inicial v0=0 en el instante t=0. La transformada de Laplace de la ecuación diferencial
que describe el comportamiento del oscilador de frecuencia propia ω0 y coeficiente de
la fuerza de rozamiento γ, es F(s)

d2xdt2+2γdxdt+ω20x=f(t)F(s)=f0a(1s2−exp(−as)s2−aexp(−as)s)1s2+2γs+ω20d2xdt2+2γdxdt+
ω02x=f(t)F(s)=f0a(1s2−exp(−as)s2−aexp(−as)s)1s2+2γs+ω02

Descomponemos la fracción de la izquierda en suma de fracciones más simples

1s(s2+2γs+ω20)=A1s+C1s+D1s2+2γs+ω20A1=1ω20C1=−1ω20D1=−2γω201s2(s2+2γs+ω20)=A2s+B2s2+C2
s+D2s2+2γs+ω20A2=−2γω40B2=1ω20C2=2γω40D2=4γ2−ω20ω401s(s2+2γs+ω02)=A1s+C1s+
D1s2+2γs+ω02A1=1ω02 C1=−1ω02
D1=−2γω021s2(s2+2γs+ω02)=A2s+B2s2+C2s+D2s2+2γs+ω02A2=−2γω04
B2=1ω02 C2=2γω04 D2=4γ2−ω02ω04

Escribimos F(s) en forma apropiada para obtener la transformada inversa de Laplace,

mirando a las tablas.

F(s)=f0a(A2s+B2s2+C2s+D2s2+2γs+ω20)−f0aexp(−as)(A2s+B2s2+C2s+D2s2+2γs+ω20)−f0exp(−as
)(A1s+C1s+D1s2+2γs+ω20)F(s)=f0a(A2s+B2s2+C2s+D2s2+2γs+ω02)−f0aexp(−as)(
A2s+B2s2+C2s+D2s2+2γs+ω02)−f0exp(−as)(A1s+C1s+D1s2+2γs+ω02)

La transformada inversa de Laplace es

x(t)=f0a(g(t)−u(t−a)g(t−a)+a⋅u(t−a)h(t−a))h(t)=A1+C1exp(−γt)cos(ωt)+(D1−C1γω)
exp(−γt)sin(ωt)g(t)=A2+B2t+C2exp(−γt)cos(ωt)+(D2−C2γω)exp(−γt)sin(ωt)ω=√
ω20−γ2 x(t)=f0a(g(t)−u(t−a)g(t−a)+a·u(t−a)h(t−a))h(t)=A1+C1exp(−γt)cos(ωt
)+(D1−C1γω)exp(−γt)sin(ωt)g(t)=A2+B2t+C2exp(−γt)cos(ωt)+(D2−C2γω)ex
p(−γt)sin(ωt)ω=ω02−γ2

Aplicamos a un oscilador de frecuencia propia ω0=4, coeficiente de la fuerza de

rozamiento proporcional a la velocidad γ=0. La fuerza dura un tiempo a=0.4 s y su valor
máximo es f0=F0/m=40.

g=0; %rozamiento
w0=4; %frecuencia propia del oscilador
F0=40; %máximo valor de la fuerza
a=0.4; %tiempo de la fuerza

A1=1/w0^2;
C1=-1/w0^2;
D1=-2*g/w0^2;

A2=-2*g/w0^4;
B2=1/w0^2;
C2=2*g/w0^4;
D2=(4*g^2-w0^2)/w0^4;
w=sqrt(w0^2-g^2);

h=@(t) A1+C1*exp(-g*t).*cos(w*t)+(D1-C1*g)*exp(-g*t).*sin(w*t)/w;
g=@(t) A2+B2*t+C2*exp(-g*t).*cos(w*t)+(D2-C2*g)*exp(-
g*t).*sin(w*t)/w;

t=linspace(0,2.4,200);
x=F0*(g(t)-g(t-a).*heaviside(t-a)-a*h(t-a).*heaviside(t-a))/a;
plot(t,x)
ylabel('x(t)')
xlabel('t')
title('Pulso F_0·(t/a)')
grid on
Obtenemos un resultado similar empleando Math Symbolic, salvo que aparece un error
cuando se describe el oscilador sin rozamiento γ=0.

>> syms w0 g s t;
>> ft=100*(1-heaviside(t-0.4))*t;
>> gs=laplace(ft);
>> Fs=gs/(s^2+2*g*s+w0^2);
>> x=ilaplace(Fs);
>> xx=subs(x,{g,w0},{0.001,4}); %rozamiento muy pequeño
>> ezplot(xx,[0 2.4])
>> grid on
>> ylabel('x(t)')
>> xlabel('t')
>> title('Pulso F_0·(t/a)')

Función rampa

Aplicamos una fuerza de la forma F0·t/a durante un tiempo a≤2, la fuerza

es F0 para t>a.
La función f(t) se define

f(t)=F(t)m=(u(t)−u(t−a))f0ta+u(t−a)f0=f0a(t−u(t−a)(t−a))f(t)=F(t)m=(u(t)−u(t−a))f
0ta+u(t−a)f0=f0a(t−u(t−a)(t−a))

Su transformada de Laplace g(s) es

g(s)=f0a(1s2−exp(−as)s2)g(s)=f0a(1s2−exp(−as)s2)

Para un oscilador de frecuencia propia ω0 y coeficiente de la fuerza de rozamiento γ, y

para las condiciones iniciales: posición inicial x0=0 y velocidad inicial v0=0 en el
instante t=0, calculamos la transformada inversa de Laplace de la función.

F(s)=f0a(1s2−exp(−as)s2)1s2+2γs+ω20F(s)=f0a(1s2−exp(−as)s2)1s2+2γs+ω02

Expresamos la siguiente fracción como suma de fracciones más simples

1s2(s2+2γs+ω20)=As+Bs2+Cs+Ds2+2γs+ω20A=−2γω40B=1ω20C=2γω40D=4γ2−ω20ω401s2(s2+2γs
+ω02)=As+Bs2+Cs+Ds2+2γs+ω02A=−2γω04 B=1ω02 C=2γω04
D=4γ2−ω02ω04

Escribimos F(s) en forma apropiada para obtener la transformada inversa de Laplace,

mirando a las tablas.

F(s)=f0ωf(As+Bs2+ω2f+Cs+D(s+γ)2+(ω20−γ2))+f0πaexp(−as)(As+Bs2+ω2f+Cs+D(s+γ)2+(ω20−γ2))x
(t)=f0ωf(g(t)+u(t−a)g(t−a))g(t)=Acos(ωft)+Bsin(ωft)+Cexp(−γt)cos(ωt)+(D−Cγω
)exp(−γt)sin(ωt)ω=√ ω20−γ2 ωf=πaF(s)=f0ωf(As+Bs2+ωf2+Cs+D(s+γ)2+(ω02
−γ2))+f0πaexp(−as)(As+Bs2+ωf2+Cs+D(s+γ)2+(ω02−γ2))x(t)=f0ωf(g(t)+u(t
−a)g(t−a))g(t)=Acos(ωft)+Bsin(ωft)+Cexp(−γt)cos(ωt)+(D−Cγω)exp(−γt)sin(
ωt)ω=ω02−γ2 ωf=πa

Representamos gráficamente x(t) para f0=10, γ=0.5, ω0=3 rad/s

g=0.5; %rozamiento
w0=3; %frecuencia propia
f0=10; %valor máximo de la fuerza
a=2; %anchura del pulso

A=-2*g/w0^4;
B=1/w0^2;
C=2*g/w0^4;
D=(4*g^2-w0^2)/w0^4;
w=sqrt(w0^2-g^2);
g=@(t) A+B*t+C*exp(-g*t).*cos(w*t)+(D-C*g)*exp(-g*t).*sin(w*t)/w;

t=linspace(0,10,400);
x=f0*(g(t)-heaviside(t-a).*g(t-a))/a;
plot(t,x)
ylabel('x(t)')
xlabel('t')
title('Pulso rampa')
 grid on

Obtenemos un resultado similar empleando Math Symbolic.

>> syms a t s g w0 f0;

>> ft=f0*(t-heaviside(t-a)*(t-a))/a;
>> gs=laplace(ft);
>> Fs=gs/(s^2+2*g*s+w0^2);
>> x=ilaplace(Fs);
>> xx=subs(x,{g w0 f0 a},{0.5 3 10 2});
>> ezplot(xx,[0 10])
>> grid on
>> ylabel('x(t)')
>> xlabel('t')
>> title('Pulso rampa')

Pulso media onda

Aplicamos una fuerza de la forma F(t)=F0·sin(ωft) durante un tiempo a=π/ωf. La
función f(t) se define

f(t)=F(t)m={f0sin(ωft)0≤t≤a0t>af(t)=F(t)m={f0sin(ωft) 0≤t≤a0 t>a

La transformada de Laplace de esta función se calcula a partir de su definición,

integrando dos veces por partes.

g(s)=∫∞0e−stf(t)⋅dt=f0∫a0e−stsin(πat)⋅dt=f0ωf1+exp(−as)ω2f+s2g(s)=∫0∞e−stf(t)·dt
=f0∫0ae−stsin(πat)·dt=f0ωf1+exp(−as)ωf2+s2

Para un oscilador de frecuencia propia ω0 y coeficiente de la fuerza de rozamiento γ, y

para las condiciones iniciales: posición inicial x0=0 y velocidad inicial v0=0 en el
instante t=0, calculamos la transformada inversa de Laplace de la función.

F(s)=f0ωf1+exp(−as)ω2f+s21s2+2γs+ω20F(s)=f0ωf1+exp(−as)ωf2+s21s2+2γs+
ω02

Expresamos la siguiente fracción como suma de fracciones más simples

1(s2+ω2f)(s2+2γs+ω20)=As+Bs2+ω2f+Cs+Ds2+2γs+ω20A=−2γ(ω20−ω2f)2+4γ2B=ω20−ω2f(ω20−ω2f)2+4γ2C
=2γ(ω20−ω2f)2+4γ2D=4γ2−ω20+ω2f(ω20−ω2f)2+4γ21(s2+ωf2)(s2+2γs+ω02)=As+Bs2+ωf2+
Cs+Ds2+2γs+ω02A=−2γ(ω02−ωf2)2+4γ2 B=ω02−ωf2(ω02−ωf2)2+4γ2
C=2γ(ω02−ωf2)2+4γ2 D=4γ2−ω02+ωf2(ω02−ωf2)2+4γ2

Escribimos F(s) en forma apropiada para obtener la transformada inversa de Laplace,

mirando a las tablas.

F(s)=f0ωf(As+Bs2+ω2f+Cs+D(s+γ)2+(ω20−γ2))+f0πaexp(−as)(As+Bs2+ω2f+Cs+D(s+γ)2+(ω20−γ2))x
(t)=f0ωf(g(t)+u(t−a)g(t−a))g(t)=Acos(ωft)+Bsin(ωft)+Cexp(−γt)cos(ωt)+(D−Cγω
)exp(−γt)sin(ωt)ω=√ ω20−γ2 ωf=πaF(s)=f0ωf(As+Bs2+ωf2+Cs+D(s+γ)2+(ω02
−γ2))+f0πaexp(−as)(As+Bs2+ωf2+Cs+D(s+γ)2+(ω02−γ2))x(t)=f0ωf(g(t)+u(t
−a)g(t−a))g(t)=Acos(ωft)+Bsin(ωft)+Cexp(−γt)cos(ωt)+(D−Cγω)exp(−γt)sin(
ωt)ω=ω02−γ2 ωf=πa

Representamos gráficamente x(t) para f0=1, γ=1, a=π, ω0=√ 3 ω0=3 rad/s

g=1; %rozamiento
f0=1; %amplitud de la fuerza
w0=sqrt(3); %frecuencia propia del oscilador
a=pi; %semiperiodo
wf=pi/a; %frecuencia de la fuerza oscilante

A=-2*g/((w0^2-wf^2)^2+4*g^2);
C=-A;
B=(w0^2-wf^2)/((w0^2-wf^2)^2+4*g^2);
D=(wf^2-w0^2+4*g^2)/((w0^2-wf^2)^2+4*g^2);
w=sqrt(w0^2-g^2);
g=@(t) A*cos(wf*t)+(B/wf)*sin(wf*t)+
C*exp(-g*t).*cos(w*t)+(D-C*g)*exp(-g*t).*sin(w*t)/w;
 t=linspace(0,10,100);
x=f0*wf*(g(t)+heaviside(t-a).*g(t-a));
plot(t,x)
ylabel('x(t)')
xlabel('t')
title('Pulso sin(t)')
grid on

El máximo valor del desplazamiento es 0.3315 en el instante 2.424

Obtenemos un resultado similar empleando Math Symbolic.

>> syms a t s g w0 f0;

>> ft=f0*(1-heaviside(t-a))*sin(pi*t/a);
>> gs=laplace(ft);
>> Fs=gs/(s^2+2*g*s+w0^2);
>> x=ilaplace(Fs);
>> xx=subs(x,{g w0 f0 a},{1 sqrt(3) 1 pi});
>> ezplot(xx,[0 10])
>> grid on
>> ylabel('x(t)')
>> xlabel('t')
>> title('Pulso seno')
Nota: La versión 2014 de MATLAB no ejecuta este código, da un mensaje de error
'Division by zero'. Sin embargo, es ejecutado sin problemas por la versión 2007

Energías Renovables Curso Interactivo de Física en Internet

Angel Franco García, Copyright © 2016
Oscilaciones libres
El estudio del oscilador armónico constituye un capítulo importante de la Física, ya que
son muchos los sistemas físicos oscilantes que se dan en la naturaleza y que han sido
producidos por el hombre.

En esta página, se estudian las oscilaciones libres tomando como modelo una partícula
de masa m unida a un muelle elástico de constante k.

Cuando una partícula se desplaza x de la posición de equilibrio, actúa sobre ella

una fuerza que es proporcional al desplazamiento x y de sentido contrario a
éste, tal como se muestra en la figura.

La ecuación del movimiento se escribe

ma=-kx

Teniendo en cuenta que la aceleración es la derivada segunda de la posición x,

expresamos la ecuación del movimiento como ecuación diferencial de segundo
orden.

d2xdt2+ω20x=0ω20=kmd2xdt2+ω02x=0 ω02=km

ω0 se denomina frecuencia propia o natural del oscilador armónico.

La ventaja de expresar las oscilaciones en términos de una ecuación diferencial

es que podemos establecer analogías entre sistemas físicos oscilantes
completamente diferentes: mecánicos, eléctricos, hidráulicos, etc.

La solución de esta ecuación diferencial es

x=Asin(ω0t)+Bcos(ω0t)v=dxdt=ω0(Acos(ω0t)−Bsin(ω0t))x=Asin(ω0t)+Bco
s(ω0t)v=dxdt=ω0(Acos(ω0t)−Bsin(ω0t))

Las constantes A y B se determinan a partir de la posición inicial x0 y la

velocidad inicial v0 en el instante t=0.

t=0{x0=Bv0=ω0Ax=v0ω0sin(ω0t)+x0cos(ω0t)t=0
{x0=Bv0=ω0Ax=v0ω0sin(ω0t)+x0cos(ω0t)
Sea un oscilador de frecuencia angular ω0=100 rad/s. Sabiendo que la partícula
parte de la posición x0=5 con velocidad inicial nula, v0=0, obtener y representar
su posición x en función del tiempo t

O bien, de forma equivalente, x=A·sin(ω0t+φ), donde A es la amplitud y φ la

fase inicial

A=
⎷ v20ω20+x20 tanφ=x0ω0v0A=v02ω02+x02 tanφ=x0ω0v0

>> syms w0 x0 v0;

>> x=dsolve('D2x+w0^2*x=0','x(0)=x0','Dx(0)=v0');
>> x=simplify(x)
x =x0*cos(t*w0) + (v0*sin(t*w0))/w0
>> xx=subs(x,{w0 x0 v0},{100 5 0})
xx =5*cos(100*t)
>> ezplot(xx,[0 0.1*pi])
>> xlabel('t')
>> ylabel('x')
>> title('oscilaciones libres')
>> grid on
Trayectoria en el espacio de las fases

El espacio de las fases nos muestra otra perspectiva del comportamiento de un

oscilador. Se representa el momento lineal (o la velocidad) v en el eje vertical,
y la posición del móvil x en el eje horizontal.

x=A·sin(ω0t+φ)
v=A·ω0·cos(ω0t+φ)

Eliminando el tiempo t en estas dos ecuaciones, obtenemos la ecuación de la

trayectoria, una elipse.

x2A2+y2(Aω0)2=1x2A2+y2(Aω0)2=1

Energía del oscilador

La característica esencial de una oscilación libre es que la amplitud se mantiene
constante y por tanto, la energía total se mantiene constante.

E=12mv2+12kx2=12mv2+12mω20x2=12m(v20+x20ω20)E=12mv2+12kx2=12
mv2+12mω02x2=12m(v02+x02ω02)

Actividades
Se introduce

 la posición inicial x0, en el control titulado Posición

 la velocidad inicial del móvil v0, en el control titulado Velocidad.
 la frecuencia angular se ha fijado en ω0 =100 rad/s.

Se pulsa en el botón Nuevo.

Posición:

Velocidad:
  Inicio
 Oscilaciones
 Amortiguadas

Oscilaciones amortiguadas
Fuerza de rozamiento proporcional a la velocidad
En esta página, estudiamos las oscilaciones amortiguadas tomando como modelo una
partícula de masa m unida a un muelle elástico de constante k que experimenta una
fuerza de rozamiento proporcional a la velocidad.

Al resolver la ecuación diferencial del movimiento se distinguirán tres casos:

armortiguadas, críticas y sobreamortiguadas.

La experiencia nos muestra que la amplitud de un cuerpo vibrante tal como un resorte o
un péndulo, decrece gradualmente hasta que se detiene.

Para explicar el amortiguamiento, supondremos que además de la fuerza elástica F=-kx,

actúa otra fuerza opuesta a la velocidad Fr=-λv, donde λ es una constante que depende
del sistema físico particular. Todo cuerpo que se mueve en el seno de un fluido viscoso
en régimen laminar experimenta una fuerza de rozamiento proporcional a la velocidad y
de sentido contrario a ésta.

La ecuación del movimiento se escribe

ma=-kx-λv

Expresamos la ecuación del movimiento en forma de ecuación diferencial, teniendo en

cuenta que la aceleración es la derivada segunda de la posición x, y la velocidad es la
derivada primera de x.
d2xdt2+2γdxdt+ω20x=0ω20=km2γ=λmd2xdt2+2γdxdt+ω02x=0 ω02=km
2γ=λm

 ω0 es la frecuencia natural o propia del oscilador

 γ es la constante de amortiguamiento, γ<ω0

La solución de la ecuación diferencial para (γ<ω0) tiene la siguiente expresión

s2+2γs+ω20=0{s1=−γ+ωis2=−γ−ωiω=√ ω20−γ2 x=exp(−γt)(Asin(ωt)+Bcos(ωt)

)v=dxdt=−γexp(−γt)(Asin(ωt)+Bcos(ωt))+ωexp(−γt)(Acos(ωt)−Bsin(ωt))s2+2γ
s+ω02=0
{s1=−γ+ωis2=−γ−ωiω=ω02−γ2x=exp(−γt)(Asin(ωt)+Bcos(ωt))v=dxdt=−γexp
(−γt)(Asin(ωt)+Bcos(ωt))+ωexp(−γt)(Acos(ωt)−Bsin(ωt))

donde A y B se determinan a partir de las condiciones iniciales: posición inicial x0 y

velocidad inicial v0 en el instante t=0

t=0{x0=Bv0=−γB+ωAx=exp(−γt)(v0+γx0ωsin(ωt)+x0cos(ωt))t=0
{x0=Bv0=−γB+ωAx=exp(−γt)(v0+γx0ωsin(ωt)+x0cos(ωt))

x=A0exp(−γt)sin(ωt+φ)A0=√ A2+B2 =√ v20+2v0γx0+x20ω20ω20−γ2 tanφ=BA=x0√ ω20−γ2 v

0+γx0x=A0exp(−γt)sin(ωt+φ)A0=A2+B2=v02+2v0γx0+x02ω02ω02−γ2tanφ=B

A=x0ω02−γ2v0+γx0

La característica esencial de las oscilaciones amortiguadas es que la amplitud de la oscilación

disminuye exponencialmente con el tiempo

Sea una oscilación amortiguada de frecuencia angular propia ω0=100 rad/s y cuya
constante de amortiguamiento γ=7.0 s-1. Sabiendo que la partícula parte de la
posición x0=5 con velocidad inicial nula, v0=0, escribir la ecuación de la oscilación
amortiguada.

La frecuencia angular de la oscilación amortiguada ω es

ω=√ 1002−72 =99.75rad/sω=1002−72=99.75 rad/s

A0=5.01
tanφ=14.25

La ecuación de la oscilación amortiguada es

x=5.01·exp(-7t)·sin(99.75t+1.5)

Creamos un script para representar la oscilación y la amplitud en función del tiempo.

Denominamos

 Coeficiente de rozamiento, γ (g),

  Frecuencia angular propia, ω0 (w0),

syms g w0 t x0 v0;
x=dsolve('D2x+2*g*Dx+w0^2*x=0','x(0)=x0','Dx(0)=v0');
xx=subs(x,{g w0 x0 v0},{7 100 5 0});
hold on
ezplot(xx,[0 0.7])
A0=sqrt((v0^2+2*v0*g*x0+x0^2*w0^2)/(w0^2-g^2));
AA=subs(A0,{g w0 x0 v0},{7 100 5 0});
xp=AA*exp(-7*t);
ezplot(xp,[0 0.7])
h=ezplot(xp,[0 0.7]);
set(h,'color','r')
h=ezplot(-xp,[0 0.7]);
set(h,'color','r')
ylim([-5 5])
grid on
title('oscilaciones amortiguadas')
hold off

La energía del oscilador amortiguado

La energía de la partícula que describe una oscilación amortiguada es la suma de la
energía cinética de la partícula y de la energía potencial del muelle elástico deformado.

E=12mv2+12kx2=12mv2+12mω20x2E=12mv2+12kx2=12mv2+12mω02x2

Introducimos las expresiones de la posición x y de la velocidad v de la partícula en

función del tiempo t.
E=12mω20A2e−2γt−12mγωA2e−2γtsin(2(ωt+φ))E=12mω02A2e−2γ t−12mγ ωA
2e−2γ tsin(2(ωt+φ))

Si la constante de amortiguamiento γ es pequeña, como hemos visto en el ejemplo del

apartado anterior ω0≈ω

E=12mω20A2e−2γt(1−γω0sin(2(ω0t+φ)))E=12mω02A2e−2γ t(1−γω0sin⁡(2(ω0t
+φ)))

La energía decrece exponencialmente con el tiempo, pero con una pequeña ondulación
debida al segundo término entre paréntesis.

Añadimos al script anterior las siguientes líneas para representar la energía del oscilador
amortiguado en función del tiempo. La energía del oscilador decrece rápidamente con el
tiempo.

%energía
v=diff(x,t);
e=0.5*v^2+0.5*100^2*x^2; %la masa m es un factor de escala m=1
ee=subs(e,{g w0 x0 v0},{7 100 5 0})
figure
ezplot(ee,[0 0.7])
grid on
title('Energía')

La energía perdida hasta el instante t a causa de la fuerza de rozamiento se calcula

mediante la integral
t∫0−(λv)v⋅dt=−2mγt∫0v2dt∫0t−(λv)v·dt=−2mγ∫0tv2dt

que será igual a la diferencia entre la energía del oscilador en el instante t y la energía
inicial del oscilador en el instante t=0

Comprobamos con MATLAB

clear
syms t g w0 x0 v0;
w=sqrt(w0^2-g^2);
x=exp(-g*t)*((v0+g*x0)*sin(w*t)/w+x0*cos(w*t)); %posición
v=diff(x,t); %velocidad
Ei=v0^2/2+w0^2*x0^2/2; %energía inicial del oscilador
Ef=v^2/2+w0^2*x^2/2; %energía final
DE=Ef-Ei; %diferencia de energías
%características del oscilador, masa m=1
E1=subs(DE,{g w0 x0 v0},{7 100 5 0});
Ep=-int(2*g*v^2,t,0,t); %trabajo de la fuerza de rozamiento
E2=subs(Ep,{g w0 x0 v0},{7 100 5 0});
subs(E1,t,0.7) %valor para el instante t=0.7 s
subs(E2,t,0.7)
Corremos el script en la ventana de comandos

ans = -1.2499e+005
ans = -1.2499e+005

Valores medios
La posición y velocidad de un oscilador libre, una partícula de masa m unida a un
muelle de constante elástica k son

x=Asin(ω0t+φ)ω20=kmv=dxdt=Aω0cos(ω0t+φ)x=Asin(ω0t+φ)
ω02=kmv=dxdt=Aω0cos(ω0t+φ)

Las energías cinética y potencial son, respectivamente

Ep=12kx2=12kA2sin2(ω0t+φ)Ek=12mv2=12mA2ω20cos2(ω0t+φ)Ep=12kx2=12kA
2sin2(ω0t+φ)Ek=12mv2=12mA2ω02cos2(ω0t+φ)

La suma de las energías cinética y potencial es constante e independiente del tiempo t

E=12kx2+12mv2=12mω20A2sin2(ω0t+φ)+12mω20A2cos2(ω0t+φ)=12mω20A2E
=12kx2+12mv2=12mω02A2sin2(ω0t+φ)+12mω02A2cos2(ω0t+φ)=12mω02
A2

Se denomina valor medio de una función periódica de periodo P a la integral

⟨f(t)⟩=1PP∫0f(t)⋅dt〈f(t)〉=1P∫0Pf(t)·dt

fplot(@(t) [sin(t), sin(t).^2],[0,2*pi])

grid on
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'})
 legend('sin(t)', 'sin(t)^2', 'location','southwest')
xlabel('t')
ylabel('x')
title('Funciones periódicas')

La función sin2(ω0t+φ) es una función periódica de periodo P=π/ω0

⟨sin2(ω0t+φ)⟩=ω0ππ/ω0∫0sin2(ω0t+φ)⋅dt=ω02ππ/ω0∫0{1−cos(2ω0t+2φ)}dt=ω02π{t−12ω
0sin(2ω0t+2φ)}π/ω00=ω02π{πω0−12ω0sin(2π+2φ)+12ω0sin(2φ)}=12〈sin2(ω0t+φ)〉=ω

0π∫0π/ω0sin2(ω0t+φ)·dt=ω02π∫0π/ω0{1−cos(2ω0t+2φ)}dt=ω02π{t−12ω0sin(
2ω0t+2φ)}0π/ω0=ω02π{πω0−12ω0sin(2π+2φ)+12ω0sin(2φ)}=12

Del mismo modo

⟨cos2(ω0t+φ)⟩=ω0ππ/ω0∫0cos2(ω0t+φ)⋅dt=ω02ππ/ω0∫0{1+cos(2ω0t+2φ)}dt=ω02π{t+12
ω0sin(2ω0t+2φ)}π/ω00=ω02π{πω0+12ω0sin(2π+2φ)−12ω0sin(2φ)}=12〈cos2(ω0t+φ)〉=
ω0π∫0π/ω0cos2(ω0t+φ)·dt=ω02π∫0π/ω0{1+cos(2ω0t+2φ)}dt=ω02π{t+12ω0si
n(2ω0t+2φ)}0π/ω0=ω02π{πω0+12ω0sin(2π+2φ)−12ω0sin(2φ)}=12

El valor medio de la energía cinética y potencial son, respectivamente

⟨12mv2⟩=E2⟨12kx2⟩=E2〈12mv2〉=E2〈12kx2〉=E2

La ecuación diferencial del oscilador amortiguado es

md2xdt2+λdxdt+kx=0md2xdt2+λdxdt+kx=0

Multiplicando por v=dx/dt

mvdvdt+kxdxdt=−λv2ddt(12mv2+12kx2)=−λv2dEdt=−λv2mvdvdt+kxdxdt=−λv2ddt(
12mv2+12kx2)=−λv2dEdt=−λv2

Tomando valores medios

d⟨E⟩dt=−2λm⟨12mv2⟩d⟨E⟩dt=−2λm12⟨E⟩d⟨E⟩⟨E⟩=−λmdtd〈E〉dt=−2λm〈12mv2〉d〈E〉dt=−
2λm12〈E〉d〈E〉〈E〉=−λmdt

Integrando entre 0 y t

⟨E⟩=⟨E⟩0exp(−2γt)〈E〉=〈E〉0exp(−2γt)

El valor medio de la energía en cada periodo, decrece exponencialmente con el tiempo.

Oscilaciones críticas (γ=ω0)

Tomamos g=w0=100;

>> x=dsolve('D2x+2*100*Dx+100^2*x=0','x(0)=5','Dx(0)=0')
x =5/exp(100*t) + (500*t)/exp(100*t)
>> ezplot(x,[0 0.02*pi])
>> title('Oscilación crítica')
>> xlabel('t')
>> ylabel('x');
>> grid on
Oscilaciones sobreamortiguadas (γ>ω0)
Tomamos g=120 y w0=100;

>> x=dsolve('D2x+2*120*Dx+100^2*x=0','x(0)=5','Dx(0)=0')
x =exp(t*(20*11^(1/2) - 120))*((15*11^(1/2))/11 + 5/2) +
(5*11^(1/2)*(11^(1/2) - 6))/(22*exp(t*(20*11^(1/2) + 120)))
>> ezplot(x,[0 0.02*pi])
>> title('oscilación sobreamortiguada')
>> xlabel('t')
>> ylabel('x');
>> grid on
Solución numérica
Expresamos la ecuación diferencial de segundo orden que describe las oscilaciones
amortiguadas en forma de dos ecuaciones diferenciales de primer orden, para
resolverlas utilizando la función ode45 de MATLAB

d2xdt2+2γdxdt+ω20x=0{dxdt=vdvdt=−2γv−ω20xd2xdt2+2γdxdt+ω02x=0{dxdt=vdv
dt=−2γv−ω02x

w0=100; %frecuencia angular propia

g=7; %rozamiento, gamma,
x0=[5,0];%posición inicial, x0, velocidad inicial, v0
tf=[0,0.3*pi]; %intervalo de tiempo

f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)];
tspan=[0 tf];
[t,x]=ode45(f,tf,x0);
%plot(t,x(:,1)) %tiempo, posición
plot(x(:,1),x(:,2)) %posición, velocidad (espacio de las fases)
xlabel('x')
ylabel('v');
title('oscilador amortiguado')
grid on
Con la sentencia anulada plot(t,x(:,1)), obtenemos una gráfica similar a la primera figura
de esta página
El espacio de las fases nos muestra otra perspectiva del comportamiento del oscilador.
Se representa el momento lineal (o la velocidad) v en el eje vertical, y la posición del
móvil x en el eje horizontal. La trayectoria es una espiral que converge hacia el origen.

Actividades
Se introduce

 La posición inicial x0, en el control titulado Posición

 La velocidad inicial del móvil v0, en el control titulado Velocidad.
 La constante de amortiguamiento γ, en el control de edición
titulado Amortiguamiento
 Se ha fijado la frecuencia angular natural del oscilador en ω0 =100 rad/s

Se pulsa el botón titulado Nuevo.

Probar con los siguientes valores de la constante de amortiguamiento γ : 5

(amortiguadas), 100 (críticas), 110 (sobreamortiguadas).

Posición: Velocidad:
Amortiguamiento:

Medida de la constante de amortiguamiento

  En el instante t, la amplitud (máximo desplamiento) es A1=Aexp(-γt)
 En el instante t+P, la amplitud es A2=Aexp(-γ(t+P))

El decrecimiento logarítmico δ nos permite calcular fácilmente la constante de

amortiguamiento γ, cuando su valor es pequeño comparado con la frecuencia natural del
oscilador ω0.

=2π(γ/ω0)√ 1−(γ/ω0)2 γ≪ω0δ≈

A1A2=exp(−γt)exp(−γ(t+P))=exp(γP)δ=lnA1A2=γP=2πγ√ ω20−γ2
2πγω0A1A2=exp(−γt)exp(−γ(t+P))=exp(γP)δ=lnA1A2=γP=2πγω02−γ2=2π(γ/
ω0)1−(γ/ω0)2γ≪ω0 δ≈2πγω0

En una experiencia hemos medido los máximos desplazamientos

 en el instante t=0, A1=5.0

 en el instante t=0.0630, A2=3.2173

δ=ln(A1/A2)=2πγ/ω0, como ω0=100 rad/s, entonces γ=7.02.

En una experiencia de laboratorio con el péndulo de Pohl se han tomado los siguientes
datos de la amplitud de una oscilación amortiguada

t (s) A

0 20

0.905 15.8

1.810 13.3

2.715 10.5

3.620 8.7

4.525 6.8

5.430 5.7
6.335 4.4

7.240 3.4

8.145 2.8

Si la amplitud disminuye exponencialmente con el tiempo, A=A0exp(-γt), o bien,

lnA=lnA0-γtque es la ecuación de una línea recta

Representamos ln(A) en función del tiempo t y realizamos un ajuste lineal de los datos.
La pendiente de la recta será la constante γ de amortiguamiento.

t=0.905*(0:9); %tiempo
A=[20,15.8,13.3,10.5,8.7,6.8,5.7,4.4,3.4,2.8]; %amplitud
plot(t,log(A),'ro','markersize',4,'markerfacecolor','r')
grid on
xlabel('t')
ylabel('log(A')
title('Oscilaciones amortiguadas')
En la ventana gráfica de MATLAB elegimos Tools/Basic fitting, activamos la
casilla linear y obtenemos para la pendiente de la recta que mejor ajusta a los datos
experimentales p1=-0.24148. Por tanto, la constante γ=0.24 s-1

Energías Renovables Curso Interactivo de Física en Internet

Angel Franco García, Copyright © 2016
  Inicio
 Oscilaciones
 Forzadas
 Tipos de fuerza

El oscilador forzado
El estado transitorio y su evolución hacia el estado estacionario
La amplitud de una oscilación amortiguada decrece con el tiempo. Al cabo de un cierto
tiempo teóricamente infinito, el oscilador se detiene en el origen. Para mantener la
oscilación es necesario aplicar una fuerza oscilante.

El oscilador forzado, o su equivalente el circuito LRC conectado a una fuente de

corriente alterna es un ejemplo que nos permite estudiar con detalle las soluciones de
una ecuación diferencial de segundo orden. Nos permite diferenciar entre estado
transitorio y estacionario. Comprender el importante fenómeno de la resonancia

Las fuerzas que actúan sobre la partícula son:

 La fuerza que ejerce el muelle, -k·x

 La fuerza de rozamiento proporcional a la velocidad λv y de sentido contrario a ésta
 La fuerza oscilante F0·cos(ωf t) de frecuencia angular ωf

La ecuación del movimiento de la partícula es

ma=-kx-λv+F0·cos(ωf t)
Expresamos la ecuación del movimiento en forma de ecuación diferencial que describe
las oscilaciones forzadas

d2xdt2+2γdxdt+ω20x=F0m cos(ωft)ω20=km2γ=λmd2xdt2+2γdxdt+ω02x=F0m
cos(ωft) ω02=km 2γ=λm

 donde ω0 es la frecuencia natural o propia del oscilador

 ωf es la frecuencia angular de la fuerza oscilante de amplitud F
 γ es la constante de amortiguamiento, γ<ω0

La solución general de la ecuación diferencial homogénea tiene la forma

x1=(Ccos(ωt)+Dsin(ωt))⋅exp(−γt)ω=√ ω20−γ2 x1=(Ccos⁡(ω t)+Dsin(ω t))·exp

⁡(−γ t) ω=ω02−γ2

Donde los coeficientes C y D se determinan a partir de las condiciones iniciales

Supondremos inicialmente que ω0≠ωf

Una solución particular de la ecuación diferencial completa tiene la forma

x2=Acos(ωf t)+Bsin(ωf t)

Obtendremos los valores de A y B haciendo que cumpla la ecuación diferencial lineal

completa

A=F(ω20−ω2f)m((ω20−ω2f)2+4γ2ω2f)B=2γωfFm((ω20−ω2f)2+4γ2ω2f)A=F(ω02−
ωf2)m((ω02−ωf2)2+4γ2ωf2) B=2γωfFm((ω02−ωf2)2+4γ2ωf2)

La solución general de la ecuación diferencial completa es la suma de la solución

general de la homogénea más la solución particular x=x1+x2.

x=(Ccos(ωt)+Dsin(ωt))exp(−γt)+Acos(ωft)+Bsin(ωft)x=(Ccos⁡(ω t)+Dsin(ω t
))exp⁡(−γ t)+Acos⁡(ωf t)+Bsin(ωf t)

El estado transitorio desaparece al cabo de cierto tiempo (teóricamente infinito) y depende de

las condiciones iniciales

El primer término, describe el estado transitorio. El segundo término, describe el estado

estacionario.

La velocidad vale

v=dxdt=⎧⎪
⎪⎨⎪
⎪⎩−Aωfsin(ωft)+Bωfcos(ωft)−γ(Ccos(ωt)+Dsin(ωt))exp(−γt)+(−ωCsin(ωt)+ω
Dcos(ωt))exp(−γt)v=dxdt={−Aωfsin(ωft)+Bωfcos(ωft)−γ(Ccos(ω t)+Dsin(ω t
))exp(−γt)+ (−ωCsin(ω t)+ωDcos(ω t))exp(−γt)
El estado estacionario, independiente de las condiciones iniciales y es el que permanece

No hay rozamiento γ=0

Cuando ω0≠ωf
La solución de la ecuación diferencial para las condiciones iniciales: en el instante t=0,
la posición es x0 y la velocidad v0, se escriben.

x=Ccos(ω0t)+Dsin(ω0t)+Fm(ω20−ω2f)cos(ωft)v=dxdt=ω0(−Csin(ω0t)+Dcos(ω0t))−F
ωfm(ω20−ω2f)sin(ωft)t=0⎧⎨⎩x0=C+Fm(ω20−ω2f)v0=ω0Dx=x0cos(ω0t)+v0ω0sin(ω0t)+F
m(ω20−ω2f)(cos(ωft)−cos(ω0t))x=Ccos(ω0t)+Dsin(ω0t)+Fm(ω02−ωf2)cos(ωft)v=
dxdt=ω0(−Csin(ω0t)+Dcos(ω0t))−Fωfm(ω02−ωf2)sin(ωft)t=0
{x0=C+Fm(ω02−ωf2)v0=ω0Dx=x0cos(ω0t)+v0ω0sin(ω0t)+Fm(ω02−ωf2)(co
s(ωft)−cos(ω0t)) .

Establecemos los siguientes valores:

 Cociente F/m=1
 Frecuencia angular propia, w0=100 rad/s y
 Frecuencia de la fuerza oscilante, wf=120 rad/s

>> syms t w0 wf F x0 v0;

>> x=dsolve('D2x+w0^2*x=F*cos(wf*t)','x(0)=x0','Dx(0)=v0');
>> x=simplify(x)
x =(w0*(F*cos(t*w0) - F*cos(t*wf)) - v0*w0^2*sin(t*w0) +
v0*wf^2*sin(t*w0))/(w0*wf^2 - w0^3) + x0*cos(t*w0)
>> xx=subs(x,{w0 wf F x0 v0},{100 120 1 0 0});
>> ezplot(xx,[0 0.2*pi])
>> title('Sin rozamiento')
>> xlabel('t')
>> ylabel('x')
>> grid on
Obtenemos pulsaciones, suma de armónicos de dos frecuencias distintas

Cuando ω0=ωf
Cuando la frecuencia de la fuerza oscilante ωf se hace igual a la frecuencia propia del
oscilador ω0, en el tercer término de la ecuación que nos da la posición x tenemos una
ideterminación del tipo 0/0.

limωf→ω0(cos(ωft)−cos(ω0t)ω20−ω2f)=tsin(ω0t)2ω0x=x0cos(ω0t)+v0ω0sin(ω0t)+Fmtsin(ω0t)2ω0li
mωf→ω0(cos(ωft)−cos(ω0t)ω02−ωf2)=tsin(ω0t)2ω0x=x0cos(ω0t)+v0ω0sin(
ω0t)+Fmtsin(ω0t)2ω0

>> syms w0 wf t;
>> limit((cos(wf*t)-cos(w0*t))/(w0^2-wf^2),wf,w0)
ans =(t*sin(t*w0))/(2*w0)
O bien, resolvemos de nuevo la ecuación diferencial con ωf=ω0. Establecemos los
siguientes valores:

 Coeficiente de rozamiento g=0;

 Frecuencia angular propia y de la fuerza oscilante, w0=100 rad/s

>> syms w0 F x0 v0;

>> x=dsolve('D2x+w0^2*x=F*cos(w0*t)','x(0)=x0','Dx(0)=v0');
>> x=simplify(x)
x =x0*cos(t*w0) + (v0*sin(t*w0))/w0 + (F*t*sin(t*w0))/(2*w0)
>> xx=subs(x,{w0 F x0 v0},{100 1 0 0})
xx =(t*sin(100*t))/200
>> ezplot(xx,[0 0.3*pi])
 >> title('Sin rozamiento')
>> xlabel('t')
>> ylabel('x')
>> grid on

La amplitud crece linealmente sin límite

Rozamiento γ<ω0
Si las condiciones iniciales son t=0, x=x0, v=v0.

C=x0−AD=1ω(v0+γC−ωfB)C=x0−AD=1ω(v0+γC−ωfB)

Las condiciones iniciales más sencillas son x=0, y dx/dt=0 en el instante t=0. La
partícula de masa m parte del origen con velocidad inicial nula.

C=−F(ω20−ω2f)m((ω20−ω2f)2+4γ2ω2f)D=−Fγ(ω20+ω2f)mω((ω20−ω2f)2+4γ2ω2f)
C=−F(ω02−ωf2)m((ω02−ωf2)2+4γ2ωf2)
D=−Fγ (ω02+ωf2)mω((ω02−ωf2)2+4γ2ωf2)

La posición de la partícula x que experimenta una oscilación forzada en función del

tiempo t es

x=Fm((ω20−ω2f)2+4γ2ω2f)⎛⎜
⎜⎝(ω20−ω2f)(cos(ωft)−exp(−γt)cos(ωt))+2γωf(sin(ωft)−ω20+ω2f2ωωfexp(−γt)sin(
ωt))⎞⎟
⎟⎠x=Fm((ω02−ωf2)2+4γ2ωf2)((ω02−ωf2)(cos⁡(ωft)−exp⁡(−γ t)cos⁡(ω t))
+2γωf(sin(ωft)−ω02+ωf22ωωfexp⁡(−γ t)sin(ω t)))

Resolvemos la ecuación diferencial, llamando

 Coeficiente de rozamiento, γ (g),

 Frecuencia angular propia, ω0 (w0),
 Frecuencia angular de la fuerza oscilante, ωf (wf)

y estableciendo el valor de la amplitud de la fuerza oscilante F/m=1.

Cuando ω0≠ωf
Establecemos los siguientes valores:

 Coeficiente de rozamiento g=7;

 Frecuencia angular propia, w0=100 rad/s y
 Frecuencia de la fuerza oscilante, wf=120 rad/s

>> syms wf w0 g F x0 v0;

>> x=dsolve('D2x+2*g*Dx+w0^2*x=F*cos(wf*t)','x(0)=0','Dx(0)=0');
>> xx=subs(x,{g w0 wf F x0 v0},{7 100 120 1 0 0});
>> ezplot(xx,[0 0.3*pi])
>> xlabel('t')
>> ylabel('x')
>> title('Con rozamiento')

Al cabo de un cierto tiempo (teóricamente infinito) el estado transitorio desaparece y la

amplitud de la oscilación forzada tiende hacia un valor constante.
Cuando ω0=ωf
Establecemos los siguientes valores:

 Coeficiente de rozamiento g=7;

 Frecuencia angular propia y de la fuerza oscilante, w0=100 rad/s

>> syms w0 g F x0 v0;

>> x=dsolve('D2x+2*g*Dx+w0^2*x=F*cos(w0*t)','x(0)=0','Dx(0)=0');
>> xx=subs(x,{g w0 F x0 v0},{7 100 1 0 0});
>> ezplot(xx,[0 0.3*pi])
>> xlabel('t')
>> ylabel('x')
>> title('Con rozamiento')

La amplitud de la oscilación crece y tiende hacia un valor límite constante

Actividades
Se introduce

 la posición inicial x0, en el control titulado Posición

 la velocidad inicial del móvil v0, en el control titulado Velocidad.
 la constante de amortiguamiento γ, en el control titulado Amortiguamiento
 la amplitud F de la fuerza oscilante, en el control titulado Amplitud
 La frecuencia angular ωf de la fuerza oscilante, en el control titulado Frecuencia
 la frecuencia angular natural del oscilador ω0=100 rad/s no se puede modificar

Se pulsa en el botón Nuevo.

Se observa la posición del móvil en función del tiempo

Ejemplos de oscilaciones forzadas

La amplitud F de la fuerza oscilante nos permite variar la escala vertical de la
representación gráfica. Si la posición x(t) crece más allá de los límites de la ventana, se
reduce el valor de la amplitud F en el control correspondiente.

La frecuencia propia del oscilador en ω0=100 rad/s, no se puede cambiar, pero se puede
variar la frecuencia de la fuerza oscilante ωf alrededor de dicha frecuencia.

Observar y describir cada una de las representaciones

Dos opciones se presentan para el estudio completo de las oscilaciones forzadas.

1. Condiciones iniciales fijadas de antemano en (x=0, v=0), el móvil se encuentra en el

origen en el instante inicial.
2. Condiciones iniciales que puede seleccionar el usuario del programa. Se trata de
comprobar que el estado transitorio depende de las condiciones iniciales, pero no el
estado estacionario (el que describe el comportamiento del oscilador, después de un
cierto tiempo, teóricamente infinito. En la práctica, un intervalo de tiempo tanto más
pequeño cuanto mayor sea la constante de amortiguamiento).

Condiciones iniciales fijadas de antemano en (x=0, v=0),

1. Cerca de la resonancia ωf=110 y 90

2. En la resonancia ωf=100
3. Observar pulsos, para ωf=80 y γ=1.

Condiciones iniciales fijadas por el usuario

 Con rozamiento (γ=10)

Ejemplos de condiciones iniciales en la resonancia (ωf =100)

o x=-5, v=0
o x=+5, v=0
o x=0, v=+500
o x=0, v=-500

Ejemplos de condiciones iniciales en las proximidades de la resonancia (ωf =90)

o x=-5, v=0
o x=+5, v=0
o x=0, v=+500
o x=0, v=-500
 Sin rozamiento (γ =0)

Ejemplos de condiciones iniciales en la resonancia (ωf =100)

o x=1.5, v=0
 o x=0, v=-50
o x=0, v=150

Ejemplos de condiciones iniciales en las proximidades de la resonancia (ωf =90)

o x=5, v=0
o x=0, v=-10

Frec. fuerza: Amplitud:

Posición: Velocidad:
Amortiguamiento:

Solución numérica
Escribimos la ecuación diferencial de segundo orden que describe las oscilaciones
forzadas en forma de dos ecuaciones diferenciales de primer orden, para resolverlas
utilizando la función ode45 de MATLAB

d2xdt2+2γdxdt+ω20x=F0mcos(ωft)⎧⎨⎩dxdt=vdvdt=−2γv−ω20x+F0mcos(ωft)d2xdt2+2
γdxdt+ω02x=F0mcos(ωft){dxdt=vdvdt=−2γv−ω02x+F0mcos(ωft)

Si a la frecuencia de la fuerza oscilante ωf le damos los valores wf=120 y wf=100

obtenemos las mismas figuras que resolviendo la ecuación diferencial de forma
analítica.

w0=100; %frecuencia angular propia

g=7; %rozamiento, gamma,
wf=100; %frecuencia de la fuerza oscilante (cambiar este valor)
F=1; %F0/m amplitud de la fuerza oscilante
%condiciones iniciales
x0=[0,-0.08]; %posición inicial, x0, velocidad inicial, v0
tf=0.3*pi; %tiempo final
f=@(t,x) [x(2);-2*g*x(2)-w0*w0*x(1)+F*cos(wf*t)];
tspan=[0 tf];
[t,x]=ode45(f,tspan,x0);
plot(t,x(:,1))
xlabel('t')
ylabel('x');
title('oscilador forzado')
En este ejemplo, ωf=100. Las condiciones iniciales son: x0=0, v0=-0.08. El estado
transitorio cambia, pero no lo hace la amplitud constante de la oscilación que se
establece al cabo de un cierto tiempo (teóricamente infinito)
Energías Renovables Curso Interactivo de Física en Internet
Angel Franco García, Copyright © 2016
  Inicio
 Oscilaciones
 Acoplados
 Acoplados y forzados

Dos osciladores acoplados

En esta página, vamos a introducir el concepto de modo normal de vibración y a aplicar
el principio de superposición estudiando un sistema formado por dos partículas unidas a
muelles elásticos.

Sea un sistema formado por dos osciladores acoplados, formado por dos partículas de
masas m1 y m2 situadas en los extremos de dos muelles de constantes elásticas k1 y k3.
El acoplamiento se efectúa uniendo las dos partículas mediante un muelle de
constante k2, tal como se puede ver en la figura.

Llamemos x1 y x2 a los desplazamientos de cada una de las partículas a partir de su

posición de equilibrio, medidos como positivos cuando están a la derecha. El muelle de
la izquierda se ha estirado x1, el de la derecha se ha comprimido x2 y el central se ha
deformado x2-x1. Las fuerzas sobre cada una de las partículas se indican en la figura.

 Sobre la partícula de la izquierda, se ejerce una fuerza hacia la izquierda k1x1 y una
fuerza hacia la derecha debida a la deformación del muelle central k2(x2-x1),
suponemos que x2 es mayor que x1.

 Sobre la partícula derecha, se ejerce una fuerza hacia la izquierda k3x2 y otra fuerza
hacia la izquierda debida a la deformación del muelle central k2(x2-x1)
El muelle central ejerce fuerzas iguales y de sentido contrario sobre cada una de las
partículas.

Aplicamos la segunda ley de Newton a cada una de las partículas, y escribimos las
ecuaciones del movimiento en forma de ecuaciones diferenciales de segundo orden

m1d2x1dt2=−k1x1+k2(x2−x1)m2d2x2dt2=−k3x2−k2(x2−x1)m1d2x1dt2=−k1x1+k2(x2
−x1)m2d2x2dt2=−k3x2−k2(x2−x1)

En forma matricial

(m100m2)d2dt2(x1x2)+((k1+k2)−k2−k2(k2+k3))(x1x2)=0Md2xdt2+Kx=0(m100m2)d
2dt2(x1x2)+((k1+k2)−k2−k2(k2+k3))(x1x2)=0Md2xdt2+Kx=0

Buscamos una solución de la forma

x1=X1sin(ωt+φ), x2=X2sin(ωt+φ)

que representa MAS de amplitud X1, X2 y frecuencia angular ω.

{(k1+k2−m1ω2)X1−k2X2=0−k2X1+(k2+k3−m2ω2)X2=0⎛⎝k1+k2m1−ω2−k2m1−k2m2k2
+k3m2−ω2⎞⎠(X1X2)=0{(k1+k2−m1ω2)X1−k2X2=0−k2X1+(k2+k3−m2ω2)X2
=0(k1+k2m1−ω2−k2m1−k2m2k2+k3m2−ω2)(X1X2)=0

Tenemos un sistema homogéneo, los cuadrados de las frecuencias de los modos

normales de vibración se calculan haciendo que el determinante de los coeficientes sea
igual a cero

∣∣
∣∣k1+k2m1−ω2−k2m1−k2m2k2+k3m2−ω2∣∣
∣∣=0|k1+k2m1−ω2−k2m1−k2m2k2+k3m2−ω2|=0

Obtenemos una ecuación de segundo grado en ω2 con dos raíces. Para cada una de las
dos frecuencias angulares ω1 y ω2 el sistema homogéneo nos proporciona una relación
entre X1 y X2 que denominaremos r1 y r2.

r1=X(1)2X(1)1=k1+k2−m1ω21k2=k2k2+k3−m2ω21r2=X(2)2X(2)1=k1+k2−m1ω22k2=k2k2+k3−m2ω22r1=X2(1
)X1(1)=k1+k2−m1ω12k2=k2k2+k3−m2ω12r2=X2(2)X1(2)=k1+k2−m1ω22k
2=k2k2+k3−m2ω22

El movimiento resultante de cada una de las partículas x1 y x2 es la combinación lineal

de los dos modos normales de vibración de frecuencias angulares ω1 y ω2

⎧⎪
⎪
⎪⎨⎪
⎪
⎪⎩x1=X(1)1sin(ω1t+φ1)+X(2)1sin(ω2t+φ2)x2=X(1)2sin(ω1t+φ1)+X(2)2sin(ω2t+φ2)
=r1X(1)1sin(ω1t+φ1)+r2X(2)1sin(ω2t+φ2){x1=X1(1)sin(ω1t+φ1)+X1(2)sin(ω2t+
φ2)x2=X2(1)sin(ω1t+φ1)+X2(2)sin(ω2t+φ2)=r1X1(1)sin(ω1t+φ1)+r2X1(2)si
n(ω2t+φ2)

A partir de las condiciones iniciales

t=0{x1=x01x2=x02⎧⎨⎩dx1dt=v01dx2dt=v02t=0 {x1=x01x2=x02
{dx1dt=v01dx2dt=v02

determinamos las constantes desconocidas: φ1,φ2,X(1)1,X(2)1φ1,φ2,X1(1),X1(2)

⎧⎪
⎪
⎪
⎪
⎪
⎪
⎪⎨⎪
⎪
⎪
⎪
⎪
⎪
⎪⎩x01=X(1)1sin(φ1)+X(2)1sin(φ2)x02=r1X(1)1sin(φ1)+r2X(2)1sin(φ2)v01=ω1X(1)1c
os(φ1)+ω2X(2)1cos(φ2)v02=ω1r1X(1)1cos(φ1)+ω2r2X(2)1cos(φ2)X(1)1sin(φ1)=r2x01−
x02r2−r1X(1)1cos(φ1)=r2v01−v02ω1(r2−r1)X(2)1sin(φ2)=r1x01−x02r1−r2X(2)1cos(φ2)=r1v01−v02ω2(r
1−r2){x01=X1(1)sin(φ1)+X1(2)sin(φ2)x02=r1X1(1)sin(φ1)+r2X1(2)sin(φ2)v01

=ω1X1(1)cos(φ1)+ω2X1(2)cos(φ2)v02=ω1r1X1(1)cos(φ1)+ω2r2X1(2)cos(φ
2)X1(1)sin(φ1)=r2x01−x02r2−r1
X1(1)cos(φ1)=r2v01−v02ω1(r2−r1)X1(2)sin(φ2)=r1x01−x02r1−r2
X1(2)cos(φ2)=r1v01−v02ω2(r1−r2)

Ejemplo
Vamos a determinar la respuesta del sistema formado por dos partículas de
masas m1=10 y m2=1 unidas por muelles de constantes k1=30, k2=5 y k3=0. Las
condiciones iniciales en el instante t=0, son las siguientes: posición inicial x01=1, x02=0,
velocidad inicial v01=0, v02=0.

% dos osciladores acoplados

k1=30; k2=5; k3=0;
m1=10; m2=1;

%cuadrado de las frecuencias de los modos normales de oscilación

a=(k1+k2)/m1+(k2+k3)/m2;
b=((k1+k2)*(k2+k3)-k2^2)/(m1*m2);
w2=roots([1,-a,b])
 r(1)=(k1+k2-m1*w2(1))/k2;
r(2)=(k1+k2-m1*w2(2))/k2;

%condiciones iniciales
X0=[1,0];
V0=[0,0];
X1_sin=(r(2)*X0(1)-X0(2))/(r(2)-r(1));
X1_cos=(r(2)*V0(1)-V0(2))/(sqrt(w2(1))*(r(2)-r(1)));
X2_sin=(r(1)*X0(1)-X0(2))/(r(1)-r(2));
X2_cos=(r(1)*V0(1)-V0(2))/(sqrt(w2(2))*(r(1)-r(2)));
X1=sqrt(X1_sin^2+X1_cos^2);
fi1=atan2(X1_sin,X1_cos);
X2=sqrt(X2_sin^2+X2_cos^2);
fi2=atan2(X2_sin,X2_cos);

%ecuación del movimiento de las partículas

t=linspace(0,20,400);
x1=X1*sin(sqrt(w2(1))*t+fi1)+X2*sin(sqrt(w2(2))*t+fi2);
x2=r(1)*X1*sin(sqrt(w2(1))*t+fi1)+r(2)*X2*sin(sqrt(w2(2))*t+fi2);

%representación gráfica
hold on
plot(t,x1,'b');
plot(t,x2,'r');
title('Dos osciladores acoplados')
ylabel('x_1,x_2')
xlabel('t')
grid on
hold off

En la ventana de comandos aparecen los cuadrados de las frecuencias de los modos

normales de vibración
 w2 =
6.0000
2.5000

Modos normales de vibración

El caso más sencillo es aquél en el que k1=k3=k y m1=m2=m.Las partículas tienen
masa m y las constantes elásticas de sus muelles son iguales a k. El muelle que acopla a
las dos partículas tiene constante k2.

Tenemos un sistema homogéneo del tipo (A-λI)X=0, donde I es la matriz unidad. Los
cuadrados de las frecuencias de los modos normales de vibración son los valores
propios λ=ω2de la matriz A.

∣∣
∣∣k+k2m−ω2−k2m−k2mk+k2m−ω2∣∣
∣∣=0|k+k2m−ω2−k2m−k2mk+k2m−ω2|=0

Calculamos mediante la función eig de MATLAB los valores propios y vectores

propios de la matriz A

(k+k2m−k2m−k2mk+k2m)(k+k2m−k2m−k2mk+k2m)

>> syms k k2 m;
>> A=[(k+k2)/m,-k2/m;-k2/m,(k+k2)/m];
>> [V,D]=eig(A)
V =
[ 1, -1]
[ 1, 1]
D =
[ k/m, 0]
[ 0, (k + 2*k2)/m]
Los valores propios λ=ω2 (cuadrado de las frecuencias angulares) están en la diagonal
de la matriz D.

ω21=kmω22=k+2k2mω12=km ω22=k+2k2m

Sus correspondientes vectores propios son los vectores columna de la matriz V.

D=(ω2100ω22)V=(X(1)1X(2)1X(1)2X(2)2)D=(ω1200ω22)
V=(X1(1)X1(2)X2(1)X2(2))

 El vector propio correspondiente al primer valor propio es

ω1→(X(1)1X(1)2)=X(1)1(V11V21)=X(1)1(11)ω1→(X1(1)X2(1))=X1(1)(V1
1V21)=X1(1)(11)

o simplemente [1;1]

 El vector propio correspondiente al segundo valor propio es

 ω2→(X(2)1X(2)2)=X(2)1(V12V22)=X(2)1(1−1)ω2→(X1(2)X2(2))=X1(2)(V
12V22)=X1(2)(1−1)

o simplemente [1;-1] ó [-1;1]

La ecuación del movimiento de cada una de las partículas es la superposición de los

modos normales de vibración

{x1=X(1)1sin(ω1t+φ1)+X(2)1sin(ω2t+φ2)x2=X(1)2sin(ω1t+φ1)+X(2)2sin(ω2t+φ2){x
1=X1(1)sin(ω1t+φ1)+X1(2)sin(ω2t+φ2)x2=X2(1)sin(ω1t+φ1)+X2(2)sin(ω2t
+φ2)

Teniendo en cuenta la relación entre amplitudes

{x1=X(1)1sin(ω1t+φ1)+X(2)1sin(ω2t+φ2)x2=X(1)1sin(ω1t+φ1)−X(2)1sin(ω2t+φ2){x
1=X1(1)sin(ω1t+φ1)+X1(2)sin(ω2t+φ2)x2=X1(1)sin(ω1t+φ1)−X1(2)sin(ω2t
+φ2)

A partir de las condiciones iniciales

t=0{x1=x01x2=x02⎧⎨⎩dx1dt=v01dx2dt=v02t=0 {x1=x01x2=x02
{dx1dt=v01dx2dt=v02

determinamos las constantes desconocidas φ1,φ2,X(1)1,X(2)1φ1,φ2,X1(1),X1(2)

⎧⎪
⎪
⎪
⎪
⎪
⎪
⎪⎨⎪
⎪
⎪
⎪
⎪
⎪
⎪⎩x01=X(1)1sin(φ1)+X(2)1sin(φ2)x02=X(1)1sin(φ1)−X(2)1sin(φ2)v01=ω1X(1)1cos(
φ1)+ω2X(2)1cos(φ2)v02=ω1X(1)1cos(φ1)−ω2X(2)1cos(φ2){x01=X1(1)sin(φ1)+X1
(2)sin(φ2)x02=X1(1)sin(φ1)−X1(2)sin(φ2)v01=ω1X1(1)cos(φ1)+ω2X1(2)cos
(φ2)v02=ω1X1(1)cos(φ1)−ω2X1(2)cos(φ2)

Que para este caso son

X(1)1sin(φ1)=x01+x022X(1)1cos(φ1)=v01+v022ω1X(2)1sin(φ2)=x01−x022X(2)1cos(φ2)=v01−v0
22ω2X1(1)sin(φ1)=x01+x022 X1(1)cos(φ1)=v01+v022ω1
X1(2)sin(φ2)=x01−x022 X1(2)cos(φ2)=v01−v022ω2

 El primer modo normal de vibración de frecuencia ω1 se obtiene cuando X(2)=0 o bien,

cuando x01=x02 y v01=v02. Los dos osciladores se mueven en fase, el muelle central no
sufre ninguna deformación y por tanto, no ejerce ninguna fuerza sobre las partículas,
las cuales se mueven como si no estuvieran acopladas.
 El segundo modo normal de vibración se obtiene cuando X(1)=0 es decir, cuando x01=-
x02 y v01=-v02. Los dos osciladores se mueven en oposición de fase.

Ejemplo
Supongamos que la constante de los muelles es k=10 y la constante del
acoplamiento k2=0.5, la masa m=1 de ambas partículas. Se desvía la primera
partícula x01=1 de la posición de equilibrio y se suelta. Vamos a determinar el
movimiento de cada partícula.

X(1)1sin(φ1)=12X(1)1cos(φ1)=0X(2)1sin(φ2)=12X(2)1cos(φ2)=0x1=12sin(ω1t+π2)+12
sin(ω2t+π2)=12cos(ω1t)+12cos(ω2t)x2=12sin(ω1t+π2)−12sin(ω2t+π2)=12cos(ω1t)−12
cos(ω2t)X1(1)sin(φ1)=12 X1(1)cos(φ1)=0 X1(2)sin(φ2)=12
X1(2)cos(φ2)=0x1=12sin(ω1t+π2)+12sin(ω2t+π2)=12cos(ω1t)+12cos(ω2t)x2
=12sin(ω1t+π2)−12sin(ω2t+π2)=12cos(ω1t)−12cos(ω2t)

% dos osciladores acoplados

k=10; k2=0.5; k=10;
m=1;

%calcula las frecuencias de los modos normales de vibración

A=[(k+k2)/m,-k2/m;-k2/m,(k+k2)/m];
[V,D]=eig(A);
%frecuencias de los modos normales de vibración
w1=sqrt(D(1,1)); w2=sqrt(D(2,2));

%calcula las amplitudes X1, X2 y

%fases fi1, fi2 a partir de las condiciones iniciales
x01=1; x02=0;
v01=0; v02=0;
X_sin1=(x01+x02)/2;
X_sin2=(x01-x02)/2;
X_cos1=(v01+v02)/(2*w1);
X_cos2=(v01-v02)/(2*w2);
X1=sqrt(X_sin1^2+X_cos1^2);
X2=sqrt(X_sin2^2+X_cos2^2);
fi1=atan2(X_sin1,X_cos1);
fi2=atan2(X_sin2,X_cos2);

%representación gráfica
t=linspace(0,40,400);
x1=X1*sin(w1*t+fi1)+X2*sin(w2*t+fi2);
x2=X1*sin(w1*t+fi1)-X2*sin(w2*t+fi2);
hold on
plot(t,x1,'b');
plot(t,x2,'r');
title('Dos osciladores acoplados')
ylabel('x_1,x_2')
xlabel('t')
 hold off

Estudio energético
La energía total del sistema es la suma de las energías cinética y potencial. Tenemos la
energía cinética de cada una de las partículas, la energía potencial elástica del muelle
izquierdo que se deforma x1, del muelle derecho que se deforma x2, y del muelle central
que se deforma x2-x1.

E=Ek+Ep=12mv21+12mv21+12kx21+12kx22+12k2(x2−x1)2E=(12mv21+12(k+k2)x21)+(12
mv22+12(k+k2)x22)−k2x1x2E=Ek+Ep=12mv12+12mv12+12kx12+12kx22+12k
2(x2−x1)2E=(12mv12+12(k+k2)x12)+(12mv22+12(k+k2)x22)−k2x1x2

Una vez agrupados los términos, el primer paréntesis depende solamente de x1 y puede
denominarse la energía del primer oscilador, el segundo término depende solamente
de x2 y puede llamarse energía del segundo oscilador. El último término, que depende
de x1 y x2 se denomina energía de acoplamiento o de interacción. Este término es el que
describe el intercambio de energía entre los dos osciladores.

Actividades
Se introduce:

 la constante elástica de los dos osciladores, en el control titulado Cte. muelles.

 la constante elástica del muelle central, en el controltitulado Cte. acoplamiento
  la masa de las partículas se ha tomado como la unidad.
 la posición inicial de la partícula de la izquierda (una cantidad menor o igual que la
unidad), en el control titulado Posición inicial de 1
 la posición inicial de la partícula de la derecha, en el control titulado Posición inicial de
2.
 las velocidades iniciales se toman como cero.

Se pulsa el botón titulado Nuevo

Modos normales de vibración

 Primer modo normal de oscilación: introducir la misma cantidad, por ejemplo, 1.0 en
los controles titulados Posición inicial 1 y Posición inicial 2.

 Segundo modo normal: introducir la misma cantidad pero con signos opuestos en
dichos controles, por ejemplo, 1.0 en Posición inicial 1, y –1.0 en Posición inicial 2.

k muelles:
k2 acoplamiento:
Posición inicial 1 Posición inicial 2

Energías Renovables Curso Interactivo de Física en Internet

Angel Franco García, Copyright © 2016