¿Por qué la solución numérica de una EDO se aleja de un equilibrio inestable?

15

Deseo simular el comportamiento de un sistema de doble péndulo. El sistema es un robot manipulador de 2 grados de libertad que no se acciona y, por lo tanto, se comportará principalmente como un doble péndulo afectado por la gravedad. La única diferencia principal con un péndulo doble es que está compuesto por dos cuerpos rígidos con propiedades de masa e inercia en sus centros de masa.

Básicamente, ode45programé bajo Matlab para resolver un sistema de EDO del siguiente tipo:

[10000M110M1200100M120M22][x˙1x˙2x˙3x˙4]=[x2V1G1x4V2G2]

donde es el ángulo del primer cuerpo con respecto a la horizontal, es la velocidad angular del primer cuerpo; es el ángulo del segundo cuerpo con respecto al primer cuerpo, y es la velocidad angular del segundo cuerpo. Todos los coeficientes se especifican en el siguiente código, en las funciones y que creé.x1x2X3X4 4rhsfMass

clear all
opts= odeset('Mass',@fMass,'MStateDependence','strong','MassSingular','no','OutputFcn',@odeplot);
sol = ode45(@(t,x) rhs(t,x),[0 5],[pi/2 0 0 0],opts);

function F=rhs(t,x)
    m=[1 1];
    l=0.5;
    a=[0.25 0.25];
    g=9.81;
    c1=cos(x(1));
    s2=sin(x(3));
    c12=cos(x(1)+x(3));
    n1=m(2)*a(2)*l;
    V1=-n1*s2*x(4)^2-2*n1*s2*x(2)*x(4);
    V2=n1*s2*x(2)^2;
    G1=m(1)*a(1)*g*c1+m(2)*g*(l*c1+a(2)*c12);
    G2=m(2)*g*a(2)*c12;

    F(1)=x(2);
    F(2)=-V1-G1;
    F(3)=x(4);
    F(4)=-V2-G2;
    F=F';     
end

function M=fMass(t,x)
    m=[1 1];
    l=0.5;
    Izz=[0.11 0.11];
    a=[0.25 0.25];
    c2=cos(x(3));
    n1=m(2)*a(2)*l;
    M11=m(1)*a(1)^2+Izz(1)+m(2)*(a(2)^2+l^2)+2*n1*c2+Izz(2);
    M12=m(2)*a(2)^2+n1*c2+Izz(2);
    M22=m(2)*a(2)^2+Izz(2);
    M=[1 0 0 0;0 M11 0 M12;0 0 1 0;0 M12 0 M22];
end

Observe cómo configuro la condición inicial de (ángulo del primer cuerpo con respecto a la horizontal) para que el sistema comience en una posición completamente vertical. De esta manera, dado que solo la gravedad está actuando, el resultado obvio es que el sistema no debería moverse en absoluto desde esa posición.X1

NOTA: en todos los gráficos a continuación, las soluciones y con respecto al tiempo.X1X3

ODE45

Cuando ejecuto la simulación durante 6 segundos ode45, obtengo la solución esperada sin ningún problema, el sistema permanece donde está y no se mueve:

ingrese la descripción de la imagen aquí

Sin embargo, cuando ejecuto la simulación durante 10 segundos, el sistema comienza a moverse sin razón:

ingrese la descripción de la imagen aquí

ODE23

Luego ejecuté la simulación ode23para ver si el problema persistía. Termino con el mismo comportamiento, solo que esta vez la divergencia comienza 1 segundo después:

ingrese la descripción de la imagen aquí

ODE15s

Luego ejecuté la simulación ode15spara ver si el problema persistía y no, el sistema parece ser estable incluso durante 100 segundos:

ingrese la descripción de la imagen aquí

Por otra parte, ode15ssolo es de primer orden y tenga en cuenta que solo hay unos pocos pasos de integración. Así que ejecuté otra simulación ode15sdurante 10 segundos pero con un MaxSteptamaño de para aumentar la precisión, y desafortunadamente, esto lleva al mismo resultado que con ambos y .0,01ode45ode23

ingrese la descripción de la imagen aquí

Normalmente, el resultado obvio de estas simulaciones sería que el sistema se mantiene en su posición inicial ya que nada lo perturba. ¿Por qué está ocurriendo esta divergencia? ¿Tiene algo que ver con el hecho de que este tipo de sistemas son de naturaleza caótica? ¿Es este un comportamiento normal para las odefunciones en Matlab?

jrojasqu
fuente
Además de las ecuaciones, creo que un esquema también ayudaría mucho a entender la pregunta.
nicoguaro
Si cree que es apropiado, puede aceptar una de las respuestas (hay un botón verde).
Ertxiem - readmitir a Monica el
No dices, pero parece que estás tramando x1y x3. (Inserte un comentario seco sobre gráficos sin leyendas o descripciones). Intente trazar los logaritmos de (los valores absolutos de) x2y x4.
Eric Towers

Respuestas:

15

Creo que Brian y Ertxiem ya han hecho los dos puntos principales: su valor inicial es un equilibrio inestable y el hecho de que sus cálculos numéricos nunca sean realmente exactos proporciona la pequeña perturbación que provocará la inestabilidad.

Para dar un poco más de detalle sobre cómo se desarrolla esto, considere su problema en forma de un problema de valor inicial general

y˙(t)=METRO-1F(t,y(t))
dondey(t)=(x1(t),x2(t),x3(t),x4(t))y

F(t,y(t))=[X2-V1-sol1X4 4-V2-sol2]

F(0 0,y0 0)=0 0y˙(0 0)=0 0F~

En su código, puede probarlo calculando

norm(rhs(0,[pi/2 0 0 0]))

lo que da 6.191e-16, por lo que casi pero no exactamente cero. ¿Cómo afecta eso a la dinámica de su sistema?

Fy0 0y~0 0

Además, en muy poco tiempo, la solución de su sistema se parece a la solución del sistema linealizado.

y˙(t)=F(0 0,y0 0)+F(0 0,y0 0)(y(t)-y0 0)=F(0 0,y0 0)(y(t)-y0 0)

FFrhsy0 0re(t): =y(t)-y0 0re

re˙(t)=F(0 0,y0 0)re(t).

No podía molestarme en calcular el jacobiano a mano, así que utilicé la diferenciación automática para obtener una buena aproximación:

J: =F(0 0,y0 0)=[0 010 00 09,810 02.45250 00 00 00 012.45250 02.45250 0]

para que tu ecuación se convierta

re˙(t)=Jre(t),re(0 0)=y~0 0-y0 0

Ahora necesitamos un paso final: podemos calcular una descomposición del valor propio del jacobiano de tal manera que

J=QreQ-1

reJQremi(t): =Q-1re(t)

mi˙(t)=remi(t),mi(0 0)=Q-1re0 0.

re

mi˙yo(t)=λyomiyo(t),miyo(0 0)=yo-th componente de Q-1re0 0

yo=1,2,3,4 4λ1=3.2485

mi1(t)=mi1(0 0)mi3.2485t.

re(0 0)=0 0mi(0 0)=Q-1re(0 0)=0 0mi1(0 0)=0 0mi1(0 0)

Daniel
fuente
16

π/ /2π/ /2

Brian Borchers
fuente
44
Si monitorea cuidadosamente las variables de estado (observando los valores impresos en notación científica), debería poder ver el movimiento inicial muy lento fuera del equilibrio.
Brian Borchers
Esto tiene sentido y, de hecho, cuando inicio el sistema en una posición vertical hacia abajo (siendo un punto de equilibrio estable), el sistema no se mueve en absoluto, al menos durante una simulación de 1000 segundos, que considero un período de tiempo muy largo .
jrojasqu
66
X1sin(0)cos(0)sin(pi/2)cos(pi/2)rhs(t,[0,0,0 0] == [0,0,0,0]
π/ /2
1
θ=0 0 10-dieciséis
4

Mire los componentes de las fuerzas calculadas en sus funciones.

π

10-dieciséis

un=1.0un=un+10-dieciséis

alephzero
fuente
4

La suposición inicial era que la posición inicial estaba en un equilibrio estable (es decir, un mínimo de la energía potencial) con energía cinética cero y el sistema comenzó a alejarse del equilibrio.
Como físicamente no puede suceder (si consideramos la mecánica clásica), me vinieron a la mente dos cosas:

  1. π/ /2-π/ /2

  2. La segunda es que quizás haya algo mal con las ecuaciones de movimiento (¿quizás un error tipográfico en alguna parte?). ¿Puedes por favor escribir las ecuaciones explícitamente? Tal vez podría trazar la aceleración angular en función de la posición inicial de cada péndulo, suponiendo que la velocidad angular sea cero para verificar si hay algo extraño.

Ertxiem - reinstalar a Monica
fuente
1
π
2
Por cierto, solo por diversión, si desea mantener el sistema en una posición vertical inestable, puede cambiar su origen de coordenadas para que el ángulo sea igual a cero apuntando hacia arriba.
Ertxiem - readmitir a Monica el
@Ertxiem otra opción es introducir una pequeña fricción en los pines que podría cometer errores numéricos.
svavil
pecado(π)
Dado que físicamente no puede suceder , dada la idea de que estamos en un equilibrio inestable, de alguna manera lo desafío. Los sistemas físicos (sin demasiada fricción) no permanecen en equilibrios inestables. En términos más generales, si simula sistemas reales, desearía evitar que se atasque accidentalmente en un equilibrio inestable (sin embargo, llegó allí); es una característica, no un error. (Hay algunas excepciones raras a esto, como el estado no infectado en inmunología, que es un equilibrio inestable que puede mantenerse).
Wrzlprmft
0

Debería buscar más sobre péndulos dobles: son lo que llamamos "sistemas caóticos". Aunque se comportan siguiendo reglas simples, a partir de condiciones iniciales ligeramente diferentes, las soluciones divergen bastante rápido. Hacer simulaciones numéricas para este tipo de sistemas no es fácil. Eche un vistazo al siguiente video para obtener más información sobre el problema.

Para el péndulo simple o doble, podría escribir una fórmula para la energía total del sistema. Suponiendo que se descuidan las fuerzas de fricción, esta energía total es preservada por el sistema analítico. Numéricamente este es un tema completamente diferente.

Antes de probar el péndulo doble, pruebe el péndulo simple. Notarás que para los métodos de Runge-Kutta de orden pequeña, la energía del sistema crecerá en las simulaciones numéricas, en lugar de permanecer constante (esto es lo que sucede en tus simulaciones: obtienes movimiento de la nada). Para evitar esto, podrían usarse métodos RK de orden superior (ode45 es de orden 4; RK de orden 8 funcionaría mejor). Hay otros métodos llamados "métodos simplécticos" que están diseñados de tal manera que las simulaciones numéricas conservan la energía. En general, debe detener la simulación tan pronto como la energía aumente significativamente en comparación con su inicialización.

Y trata de entender el péndulo simple antes de pasar al doble.

Beni Bogosel
fuente
2
Sin embargo, no se trata de que el sistema sea caótico. También puede tener un equilibrio inestable en sistemas no caóticos, por ejemplo, el péndulo simple está "sobre su cabeza", y exhibirá el mismo comportamiento descrito en la pregunta.
Daniel
1
Tampoco es cierto que la energía aumente para RKM de orden pequeño: Euler implícito es de primer orden y muestra exactamente el comportamiento opuesto.
Daniel
@BeniBogosel Usted menciona métodos simplécticos que me llamaron la atención porque obviamente, en mi ejemplo, la energía no se conserva. Sin embargo, ¿podría indicar un método simpléctico específico que podría implementarse aquí?
jrojasqu
@jrojasqu, ¿por qué dice que la energía no se conserva en su sistema?
Ertxiem - readmitir a Monica el
ode45π