Aproximación del sistema de resorte-carrito-péndulo

Question

Estoy trabajando en una simulación basada en JavaScript web de un sistema de resorte-carrito-péndulo para aplicar algunas matemáticas que he aprendido recientemente en un proyecto pequeño y ordenado. Encontré una imagen de muestra en internet de lo que estoy tratando de modelar:

Entiendo cómo derivar las ecuaciones de movimiento encontrando el Lagrangiano y luego derivando las ecuaciones de movimiento de Euler-Lagrange a partir de ahí:

Dado que la masa del carrito es $M$ y la masa de la bola en el péndulo es $m$, la longitud del péndulo $l$, la distancia del carrito desde el equilibrio del resorte (con constante $k$) $x$ y el ángulo del péndulo $\theta$ tenemos las siguientes ecuaciones de movimiento:

$$x'' cos(\theta) + l\theta'' - gsin(\theta) = 0$$ $$(m+M)x''+ml[\theta''cos(\theta) - (\theta')^2sin(\theta)] = -kx$$ (Una derivación se puede encontrar aquí, aunque este no es el punto de la pregunta).

He aislado las ecuaciones para $x''$ y $\theta''$ de la siguiente manera:

$$x'' = \frac{-kx - ml[-sin(\theta)(\theta')^2-\frac{g}{l} cos(\theta)sin(\theta)]}{m+M - mcos^2(\theta)}$$

$$\theta'' = \frac{-gsin(\theta)-cos(\theta)x''}{l}$$

Luego utilicé el método de Euler usando condiciones previas/iniciales para resolver $x''$ y luego para $\theta''$ (ya que depende de $x''$. Luego resolví $x' = x''t + x'$, $x = x't + x$ (igual para $\theta$), donde $t$ es el tamaño del paso.

Este enfoque parece funcionar bien para condiciones iniciales donde $\theta$ o $\theta'$ son pequeños (aún tengo que probar jugando con $x$). Sin embargo, cuando establezco los valores de $\theta$ en valores más grandes, entonces mi resorte se vuelve loco y va por toda la pantalla.

Nunca he configurado físicamente un sistema de resorte-carrito-péndulo, por lo que tal vez mi simulación sea correcta (aunque tengo mis dudas), pero me preguntaba qué otros métodos puedo usar para una aproximación numérica y más precisa. Siempre podría establecer el tamaño del paso para la aproximación de Euler en algo más pequeño, pero eso afectaría la velocidad y funcionalidad de mi simulación.

He escuchado hablar de cosas como Runge-Kutta pero no estoy seguro de cómo aplicarlo a un sistema de EDO de segundo orden. ¿Alguien tiene algún consejo, lectura o ejemplos que pueda consultar para aprender una aproximación más precisa?

Gracias por su tiempo.

Saludos,

Answer 1

Siempre puedes transformar un sistema de orden superior en un sistema de primer orden donde la suma de los órdenes del primero da la dimensión del segundo.

Aquí podrías definir un vector de 4 dimensiones $y$ a través de la asignación $y=(x,θ,x',θ')$ de modo que $y_1'=y_3$, $y_2'=y_4$, $x_3'=$ la ecuación para $x''$ y $y_4'=$ la ecuación para $θ''$.

En python, como el BASIC moderno, uno podría implementar esto de la siguiente manera

def prime(t,y):
    x,tha,dx,dtha = y
    stha, ctha = sin(tha), cos(tha)

    d2x = -( k*x - m*stha*(l*dtha**2 + g*ctha) ) / ( M + m*stha**2)
    d2tha = -(g*stha + d2x*ctha)/l  

    return [ dx, dtha, d2x, d2tha ]

y luego pasarlo a un solucionador de EDO por alguna variante de

sol = odesolver(prime, times, y0)

Los solucionadores disponibles son

from scipy.integrate import odeint
import numpy as np

times = np.arange(t0,tf,dt)
y0 = [ x0, th0, dx0, dth0 ]

# odeint tiene los argumentos de la función ODE intercambiados
sol = odeint(lambda y,t: prime(t,y), y0, times)

O puedes construir el tuyo

def oderk4(prime, times, y):
    f = lambda t,y: np.array(prime(t,y)) # para obtener un objeto vectorial
    sol = np.zeros_like([y]*np.size(times))
    sol[0,:] = y[:]
    for i,t in enumerate(times):
        dt = times[i+1]-t
        k1 = dt * f( t       , y        )
        k2 = dt * f( t+0.5*dt, y+0.5*k1 )
        k3 = dt * f( t+0.5*dt, y+0.5*k2 )
        k4 = dt * f( t+    dt, y+    k3 )
        sol[i+1,:] = y = y + (k1+2*(k2+k3)+k4)/6

    return sol

y luego llamarlo e imprimir la solución como

sol = oderk4(prime, times, y0)

for k,t in enumerate(times):
    print "%3d : t=%7.3f x=%10.6f theta=%10.6f" % (k, t, sol[k,0], sol[k,1])

O graficarlo

import matplotlib.pyplot as plt

x=sol[:,0]; tha=sol[:,1];
ballx = x+l*np.sin(tha); bally = -l*np.cos(tha);
plt.plot(ballx, bally)
plt.show()

donde con valores iniciales y constantes

m=1; M=3; l=10; k=0.5; g=9
x0 = -3.0; dx0 = 0.0; th0 = -1.0; dth0 = 0.0;
t0=0; tf = 150.0; dt = 0.1

Obtengo este gráfico artístico:

Answer 2

En JavaScript se puede implementar una variante universal del paso de Euler o RK4 utilizando la operación de array axpy introducida por BLAS/LAPACK

function axpy(a,x,y) { 
// devuelve a*x+y para arrays x e y de la misma longitud
    var k = y.length >>> 0;
    var res = new Array(k);
    while(k-->0) { res[k] = y[k] + a*x[k]; }
    return res;
}

Esto asume ciegamente que los arrays de entrada tienen la misma longitud.

function EulerStep(t,y,h) {
    var k = odefunc(t,y); 
    // ynext = y+h*k
    return axpy(h,k,y);
}

function RK4Step(t,y,h) {
    var k0 = odefunc(t      ,               y );
    var k1 = odefunc(t+0.5*h, axpy(0.5*h,k0,y));
    var k2 = odefunc(t+0.5*h, axpy(0.5*h,k1,y));
    var k3 = odefunc(t+    h, axpy(    h,k2,y));
    // ynext = y+h/6*(k0+2*k1+2*k2+k3);
    return axpy(h/6,axpy(1,k0,axpy(2,k1,axpy(2,k2,k3))),y);
}

La función ODE sería correspondientemente

function odefunc(t,y) {  
    var x = y[0], tha = y[1], dx = y[2], dtha = y[3];  
    var stha = Math.sin(tha),  ctha = Math.cos(tha);

    var d2x = -( k*x - m*stha*(l*dtha**2 + g*ctha) ) / ( M + m*stha**2)  
    var d2tha = -(g*stha + d2x*ctha)/l  

    return [ dx, dtha, d2x, d2tha ]
}

Ahora se puede hacer un paso de integración por fotograma, o múltiples pasos temporales con un paso de tiempo reducido en consecuencia.

La variante más complicada es usar un paso de tiempo óptimo para un nivel de error requerido e interpolar los estados para los puntos de tiempo de los fotogramas, lo que podría significar que no hay paso de integración durante varios fotogramas. Básicamente esto es lo que hacen los solucionadores dopri5/rk45/ode45 o lsoda, utilizan pasos de tiempo computados internamente e interpolan los valores para la lista de tiempos dada.