Oscilador no lineal

Sea una partícula de masa m unida a dos muelles elásticos iguales de constante k, supondremos que están en el plano horizontal y la partícula desliza en este plano sin rozamiento alguno. l0 es la longitud de los muelles sin deformar, cuando la separación entre los extremos de los muelles es 2l0.

Separamos los extremos de los muelles 2(l0+δ) y desplazamos la partícula al punto de coordendas (x,y)

La energía potencial elástica del sistema formado por los dos muelles iguales es

U(x,y)= 1 2 k ( l 1 l 0 ) 2 + 1 2 k ( l 2 l 0 ) 2 l 1 2 = ( l 0 +δy ) 2 + x 2 l 2 2 = ( l 0 +δ+y ) 2 + x 2

Las ecuaciones del movimiento de la partícula son

m d 2 x d t 2 = U x m d 2 y d t 2 = U y m d 2 x d t 2 =2kx+kx( l 0 ( l 0 +δy ) 2 + x 2 + l 0 ( l 0 +δ+y ) 2 + x 2 ) m d 2 y d t 2 =2ky+k l 0 ( l 0 +δ+y ( l 0 +δ+y ) 2 + x 2 l 0 +δy ( l 0 +δy ) 2 + x 2 )

En principio, δ podría ser positivo o negativo. En el caso de que fuese negativo, la partícula estaría en un estado inestable en el origen, por tanto, no vamos a tratar este caso. Consideraremos solamente el caso δ≥0.

Dada la posición inicial y la velocidad inicial de la partícula, se resuelve el sistema de dos ecuaciones diferenciales por procedimientos numéricos.

Casos particulares:

En esta página, estudiamos lasoscilaciones de la partícula a lo largo del eje X, que son no armónicas, incluso cuando el desplazamiento x es pequeño

Para encontrar la ecuación del movimiento de la partícula a lo largo del eje X, empezaremos por el caso más simple, δ=0, la ecuación del movimiento se reduce a

d 2 x d t 2 k m l 0 2 x 3

Oscilaciones a lo largo del eje X, δ=0

Sea un bloque de masa m unido a dos muelles elásticos iguales de constante k. La longitud de los muelles si deformar es l0, cuando el sistema está en la posición de equilibrio estable x=0.

Cuando se separa el bloque una distancia x de la posición de equilibrio y se suelta, empieza a oscilar con un periodo P.

Ecuación del movimiento

La fuerza F que ejerce cada uno de los muelles sobre el bloque es

F=k( x 2 + l 0 2 l 0 )

La resultante es

R=2Fcosθ=2F x x 2 + l 0 2 =2k( 1 l 0 x 2 + l 0 2 )x

La ecuación del movimiento es

m d 2 x d t 2 =2k( 1 l 0 x 2 + l 0 2 )x

Se transforma en

d 2 x d t 2 + 2k m x{ 1 ( 1+ ( x l 0 ) 2 ) 1/2 }

Esta ecuación, se resuelve aplicando procedimientos numéricos con las condiciones iniciales siguientes: en el instante inicial t=0, x=A, y v=dx/dt=0, siendo A la amplitud de la oscilación.

Utilizamos el procedimiento ode45 de MATLAB para resolver numéricamente la ecuación diferencial de segundo orden, para los valores de k/m=60, l0=1 y una amplitud A=0.3

k=60; %cociente k/m
l0=1; %longitud del muelle sin deformar
f=@(t,x) [x(2);-2*k*x(1)*(1-1/sqrt(1+(x(1)/l0)^2))]; 
tspan=[0 25];
[t,x]=ode45(f,[0 25],[0.3,0]);
plot(t,x(:,1),'r')
grid on
xlabel('t')
ylabel('x');
title('Dos muelles')

Solución analítica

Cuando x<<l0 la ecuación diferencial se aproxima a

d 2 x d t 2 + k m l 0 2 x 3 =0

Las condiciones iniciales son: en el instante t=0, la partícula parte de la posición x=A, con velocidad nula, v=dx/dt=0

d 2 x d t 2 = dv dt = dv dx · dx dt =v dv dx v dv dx = k m l 0 2 x 3

Integramos esta ecuación diferencial

0 v vdv = k m l 0 2 A x x 3 dx v= k 2m l 0 2 ( A 4 x 4 )

es el módulo de la velocidad

Como la partícula parte de x=+A y se dirige hacia el origen la velocidad es negativa en el primer cuarto de periodo

v= k 2m l 0 2 ( A 4 x 4 )

Para obtener la posición x del oscilador en función del tiempo t tenemos que integrar nuevamente

dx dt = k 2m l 0 2 ( A 4 x 4 ) t= 2m l 0 2 k A x dx A 4 x 4

Hacemos el cambio de variable

x=A 1 z 2

La integral se transforma

t= 1 A 2m l 0 2 k 0 z dz 1 z 2 2 z 2 A l 0 k m t= 0 z dz 1 z 2 1 1 2 z 2

Se hace el cambio de variable z=sinφ, obtenemos la integral elíptica F(θ, 1 2 )

A l 0 k m t= 0 θ dφ 1 1 2 sin 2 φ A l 0 k m t=u snu=sinθ=z

Las funciones sn y cn están relacionadas

sn2u+cn2u=1

Deshacemos el cambio de variable de z a x.

x=A 1 z 2 =A 1 sn 2 u =Acnu x=A·cn( A l 0 k m t, 1 2 )

Tomando l0=1, k/m=60, representamos el desplazamiento x del oscilador en función del tiempo t, para dos valores de la amplitud A=0.15 y A=0.3

k_m=sqrt(60);
A=0.15;
t=0:0.1:25;
[sn,cn,dn]=ellipj(A*k_m*t,1/2);
hold on
plot(t,cn*A,'b')
A=0.3;
[sn,cn,dn]=ellipj(A*k_m*t, 1/2);
plot(t,cn*A,'r')
hold off
xlabel('t')
ylabel('x')
legend ('0.15','0.3')
title('Dos muelles')
grid on

Periodo del movimiento

El tiempo que tarda en describir un cuarto de periodo P es

P 4 = 2m l 0 2 k A 0 dx A 4 x 4

Haciendo el cambio de variable x=Acosφ

P 4 = l 0 A m k 0 π/2 dφ 1 1 2 sin 2 φ P 4 = l 0 A m k K( 1 2 )

El último término es la integral elíptica completa de primera especie. Su valor aproximado es 1.8541.

>> k2=1/2;
>> ellipke(k2)
ans =    1.8541

P 7.4164 l 0 A m k

En la figura, se muestra cómo el periodo P depende de la amplitud A

Oscilaciones a lo largo del eje X, δ≠0

La ecuación del movimiento de la particula a lo largo del eje X, con x<<(l0+δ), es

d 2 x d t 2 =AxB x 3 { A=2 k m ( δ l 0 +δ ) B= k m l 0 ( l 0 +δ ) 3

Las condiciones iniciales son: en el instante t=0, la partícula parte de la posición x=0 (origen), con velocidad, v=dx/dt=v0

d 2 x d t 2 = dv dt = dv dx · dx dt =v dv dx v dv dx =AxB x 3

Integramos esta ecuación diferencial

v 0 v vdv = 0 x ( Az+B z 3 )dx ( dx dt ) 2 = 1 2 B( 2 v 0 2 B 2A B x 2 x 4 ) B 2 t= 0 x dz 2 v 0 2 B 2A B z 2 z 4 B 2 t= 0 x dz ( a 2 + z 2 )( b 2 z 2 ) { b 2 = A 2 B 2 + 2 v 0 2 B A B a 2 = A 2 B 2 + 2 v 0 2 B + A B

En el libro titulado 'Table of Integrals, Series, and Products' buscamos la solución a integrales de este tipo y la encontramos en el apartado 3.152, n° 3 de la página 279

0 x dz ( z 2 + a 2 )( b 2 z 2 ) = 1 a 2 + b 2 F( γ,r ){ γ=arcsin( x b a 2 + b 2 a 2 + x 2 ) r= b a 2 + b 2 ,bx>0

Despejamos la posición x en función del tiempo t, utilizando la función elíptica sn de Jacobi

B a 2 + b 2 2 t=F( γ,r ) sinγ=sn( B a 2 + b 2 2 t,r ) x b a 2 + b 2 a 2 + x 2 =sn( αt,r ),α= B a 2 + b 2 2 x= ab·sn( αt,r ) a 2 + b 2 cn 2 ( αt,r )

Se ha utilizado la relación

sn2u+cn2u=1

Periodo

En un Movimiento Armónico Simple, x=Asin(ωt+φ), es una función periódica de periodo P, tal que cuando se incrementa el tiempo t en P, la fase (argumento del seno) se incrementa en 2π

ω(t+P)+φ=ωt+φ+2π, P=2π/ω

Del mismo modo sn(αt,r) y cn(αt,r) son funciones periódicas que cuando el tiempo se incrementa en P, el argumento de se incrementa en 4K(r)

α(t+P)=αt+4K(r) P= 4K(r) α = 8K(r) B a 2 + b 2

Representamos la posición x en función del tiempo t durante dos periodos, para δ=0.05, el cociente k/m=60, la longitud del muelle sin deformar l0=1. La posición inicial x0=0, y la velocidad inicial es v0=0.5.

Comparamos la solución analítica y numérica, vemos que coinciden

k=60; %cociente k/m
l0=1; %longitud del muelle sin deformar
delta=0.05; %deformación inicial
v0=0.5; %velocidad inicial

A=2*k*delta/(l0+delta);
B=k*l0/(l0+delta)^3;
b2=-A/B+sqrt(A^2/B^2+2*v0^2/B);
a2=A/B+sqrt(A^2/B^2+2*v0^2/B);
r=sqrt(b2)/sqrt(a2+b2);
alfa=sqrt(B*(a2+b2)/2);
P=4*ellipke(r^2)/alfa; %periodo

t=linspace(0,2*P,200);
[sn,cn,dn]=ellipj(alfa*t,r^2);
x=sqrt(a2*b2)*sn./sqrt(a2+b2*cn.^2);
%procedimiento numérico
f=@(t,x) [x(2);-A*x(1)-B*x(1)^3]; 
[tt,xx]=ode45(f,[0,2*P],[0,v0]);
hold on
plot(t,x); %analítico
plot(tt,xx(:,1)) %numérico
hold off

grid on
xlabel('t')
ylabel('x');
title('Oscilador anarmónico')

Velocidad

Dada la posición x(t) respecto del tiempo, derivamos para calcular la velocidad, dx/dt. El lector puede consultar las fórmulas de la derivada de las funciones elípticas de Jacobi, sn(x) y cn(x)

v= dx dt =ab α·cn·dn a 2 + b 2 cn 2 + b 2 α·cn·sn·dn a 2 + b 2 cn 2 a 2 + b 2 cn 2 = αab· a 2 + b 2 ( a 2 + b 2 cn 2 ( αt,r ) ) 3/2 cn( αt,r )dn( αt,r )

Representamos la posición v en función del tiempo t durante dos periodos, para δ=0.05, el cociente k/m=60, la longitud del muelle sin deformar l0=1. La posición inicial x0=0, y la velocidad inicial es v0=0.5.

Comparamos la solución analítica y numérica, vemos que coinciden

k=60; %cociente k/m
l0=1; %longitud del muelle sin deformar
delta=0.05; %deformación inicial
v0=0.5; %velocidad inicial

A=2*k*delta/(l0+delta);
B=k*l0/(l0+delta)^3;
b2=-A/B+sqrt(A^2/B^2+2*v0^2/B);
a2=A/B+sqrt(A^2/B^2+2*v0^2/B);
r=sqrt(b2)/sqrt(a2+b2);
alfa=sqrt(B*(a2+b2)/2);
P=4*ellipke(r^2)/alfa;

t=linspace(0,2*P,200);
[sn,cn,dn]=ellipj(alfa*t,r^2);
v=alfa*sqrt(a2*b2)*(a2+b2)*cn.*dn./(a2+b2*cn.^2).^(3/2); %velocidad
%procedimiento numérico
f=@(t,x) [x(2);-A*x(1)-B*x(1)^3]; 
[tt,xx]=ode45(f,[0,2*P],[0,v0]);
hold on
plot(t,v); %analítico
plot(tt,xx(:,2)) %numérico
hold off

grid on
xlabel('t')
ylabel('v');
title('Oscilador anarmónico')  

Referencias

Detcheva V., Spassov V., A simple nonlinear oscillator: analytical and numerical solution. Phys. Educ. 28 (1993) pp. 39-42

I.S. Gradshteyn, I.M. Ryzhik. Table of Integrals, Series, and Products. Seventh Edition. Elsevier, 2007

John Thomchick, J. P. McKelvey. Anharmonic vibration of an 'ideal' Hooke's law oscillator. Am. J. Phys. 46 (1) jan. 1978, pp. 40-45