Dlaczego numeryczne rozwiązanie ODE odchodzi od niestabilnej równowagi?

Chciałbym zasymulować zachowanie systemu podwójnie wahadłowego. System jest manipulatorem robota o 2 stopniach swobody, który nie jest uruchamiany i dlatego będzie zachowywał się głównie jak podwójne wahadło pod wpływem grawitacji. Jedyną główną różnicą w przypadku podwójnego wahadła jest to, że składa się on z dwóch sztywnych ciał o właściwościach masy i bezwładności w ich środkach masy.

Zasadniczo programowałem ode45pod Matlabem, aby rozwiązać system ODE następujących typów:

$$\left[ \begin{array}{cccc} 1 & 0 & 0 & 0\\ 0 & M_{11} & 0 & M_{12}\\ 0 & 0 & 1 & 0\\ 0 & M_{12} & 0 & M_{22} \end{array} \right] \left[ \begin{array}{c} \dot{x}_1\\ \dot{x}_2\\ \dot{x}_3\\ \dot{x}_4 \end{array} \right]= \left[ \begin{array}{c} x_2\\ -V_1-G_1\\ x_4\\ -V_2-G_2 \end{array} \right]$$

gdzie $x_1$ jest kątem pierwszego ciała w stosunku do poziomu, $x_2$ jest prędkością kątową pierwszego ciała; $x_3$ to kąt drugiego ciała w stosunku do pierwszego ciała, a $x_4$ to prędkość kątowa drugiego ciała. Wszystkie te współczynniki są określone w kodeksie, przy użyciu rhsi fMassfunkcji stworzyłem.

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

Zauważ, jak ustawiłem warunek początkowy $x_1$ (kąt pierwszego korpusu w stosunku do poziomu), aby system zaczął działać w pozycji całkowicie pionowej. W ten sposób, ponieważ działa tylko grawitacja, oczywistym rezultatem jest to, że system w ogóle nie powinien się poruszać z tej pozycji.

UWAGA: we wszystkich poniższych rysunkach przedstawiłem rozwiązania $x_1$ i $x_3$ w odniesieniu do czasu.

ODE45

Kiedy uruchamiam symulację przez 6 sekund ode45, otrzymuję oczekiwane rozwiązanie bez żadnych problemów, system pozostaje tam, gdzie jest i nie porusza się:

Jednak gdy uruchamiam symulację przez 10 sekund, system zaczyna się poruszać nieracjonalnie:

ODE23

Następnie przeprowadziłem symulację, ode23aby sprawdzić, czy problem nadal występuje. Skończyłem z tym samym zachowaniem, tylko tym razem dywergencja zaczyna się 1 sekundę później:

ODE15s

Następnie uruchomiłem symulację, ode15saby sprawdzić, czy problem nadal występuje i nie, system wydaje się być stabilny nawet przez 100 sekund:

Z drugiej strony, ode15sjest tylko pierwsze zamówienie i zauważ, że jest tylko kilka kroków integrujących. Uruchomiłem więc kolejną symulację w ode15sciągu 10 sekund, ale MaxSteprozmiar $0.01$ aby zwiększyć precyzję, i niestety, prowadzi to do tego samego wyniku, co w przypadku obu ode45i ode23.

Zwykle oczywistym rezultatem tych symulacji byłby fakt, że system pozostaje w początkowej pozycji, ponieważ nic go nie zakłóca. Dlaczego występuje ta rozbieżność? Czy ma to coś wspólnego z faktem, że tego rodzaju systemy mają charakter chaotyczny? Czy jest to normalne zachowanie odefunkcji w Matlabie?

Myślę, że dwie główne kwestie zostały już przedstawione przez Briana i Ertxiema: twoja początkowa wartość jest niestabilną równowagą, a fakt, że twoje obliczenia liczbowe nigdy nie są naprawdę dokładne, zapewnia małe zaburzenie, które wywoła niestabilność.

Aby podać nieco więcej szczegółów, jak to się dzieje, rozważ swój problem w postaci ogólnego problemu z wartością początkową

$$\begin{equation} \dot{\mathbf{y}}(t) = \mathbf{M}^{-1} \mathbf{f}(t, \mathbf{y}(t)) \end{equation}$$ $\mathbf{y}(t) = (x_1(t), x_2(t), x_3(t), x_4(t))$

$$\begin{equation} \mathbf{f}(t, \mathbf{y}(t)) = \begin{bmatrix} x_2 \\ -V_1 - G_1 \\ x_4 \\ -V_2 - G_2 \end{bmatrix} \end{equation}$$

$\mathbf{f}(0, \mathbf{y}_0) = 0$$\dot{\mathbf{y}}(0) = 0$$\tilde{\mathbf{f}}$

W kodzie możesz to sprawdzić, wykonując obliczenia

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

co daje 6,191e-16 - czyli prawie, ale nie dokładnie zero. Jak to wpływa na dynamikę twojego systemu?

$\mathbf{f}$$\mathbf{y}_0$$\tilde{\mathbf{y}}_0$

Ponadto, w bardzo krótkim czasie, rozwiązanie twojego systemu wygląda jak rozwiązanie systemu zlinearyzowanego

$$\begin{equation} \dot{\mathbf{y}}(t) = \mathbf{f}(0, \mathbf{y}_0) + \mathbf{f}'(0, \mathbf{y}_0) \left( \mathbf{y}(t) - \mathbf{y}_0 \right) = \mathbf{f}'(0, \mathbf{y}_0) \left( \mathbf{y}(t) - \mathbf{y}_0 \right) \end{equation}$$

$\mathbf{f}'$$\mathbf{f}$rhs$\mathbf{y}_0$$\mathbf{d}(t) := \mathbf{y}(t) - \mathbf{y}_0$$\mathbf{d}$

$$\begin{equation} \dot{\mathbf{d}}(t) = \mathbf{f}'(0, \mathbf{y}_0) \mathbf{d}(t). \end{equation}$$

Nie mogłem sobie pozwolić na ręczne obliczenie jakobianu, więc użyłem automatycznego różnicowania, aby uzyskać dobre przybliżenie:

$$\begin{equation} \mathbf{J} := \mathbf{f}'(0, \mathbf{y}_0) = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 9.81 & 0 & 2.4525 & 0 \\ 0 & 0 & 0 & 1 \\ 2.4525 & 0 & 2.4525 & 0 \end{bmatrix} \end{equation}$$

tak aby twoje równanie stało się

$$\begin{equation} \dot{\mathbf{d}}(t) = \mathbf{J} \mathbf{d}(t), \mathbf{d}(0) = \mathbf{\tilde{y}}_0 - \mathbf{y}_0 \end{equation}$$

Teraz potrzebujemy ostatniego kroku: możemy obliczyć rozkład wartości własnej jakobianu w taki sposób, aby

$$\begin{equation} \mathbf{J} = \mathbf{Q} \mathbf{D} \mathbf{Q}^{-1} \end{equation}$$

$\mathbf{D}$$\mathbf{J}$$\mathbf{Q}$$\mathbf{d}$$\mathbf{e}(t) := \mathbf{Q}^{-1} \mathbf{d}(t)$

$$\begin{equation} \dot{\mathbf{e}}(t) = \mathbf{D} \mathbf{e}(t), \mathbf{e}(0) = \mathbf{Q}^{-1} \mathbf{d}_0. \end{equation}$$

$\mathbf{D}$

$$\begin{equation} \dot{e}_i(t) = \lambda_i e_i(t), e_i(0) = i-\textrm{th component of}~\mathbf{Q}^{-1} \mathbf{d}_0 \end{equation}$$

$i=1,2,3,4$$\lambda_1 = 3.2485$

$$\begin{equation} e_1(t) = e_1(0) e^{3.2485 t}. \end{equation}$$

$\mathbf{d}(0) = 0$$\mathbf{e}(0) = \mathbf{Q}^{-1} \mathbf{d}(0) = 0$$e_1(0) = 0$$e_1(0)$

$\pi/2$$\pi/2$

Spójrz na składowe sił obliczonych w twoich funkcjach.

$\pi$

$10^{-16}$

$a = 1.0$$a = a + 10^{-16}$

Początkowe założenie było takie, że pozycja początkowa znajdowała się w stabilnej równowadze (tj. Minimum energii potencjalnej) przy zerowej energii kinetycznej, a układ zaczął się oddalać od równowagi.
Ponieważ fizycznie tak się nie stanie (jeśli weźmiemy pod uwagę mechanikę klasyczną), przyszło mi do głowy dwie rzeczy:

1. $\pi/2$$-\pi/2$

2. Po drugie, być może coś jest nie tak z równaniami ruchu (może literówka gdzieś?). Czy możesz napisać równania wprost? Być może mógłbyś wykreślić przyspieszenie kątowe jako funkcję początkowej pozycji każdego wahadła, zakładając zerową prędkość kątową, aby sprawdzić, czy jest coś dziwnego.

Powinieneś poszukać więcej o podwójnych wahadłach: są to tak zwane „systemy chaotyczne”. Mimo że zachowują się według prostych zasad, zaczynając od nieco innych warunków początkowych, rozwiązania różnią się dość szybko. Wykonywanie symulacji numerycznych dla tego rodzaju systemów nie jest łatwe. Obejrzyj poniższy film, aby uzyskać lepszy wgląd w problem.

Dla prostego lub podwójnego wahadła możesz napisać wzór na całkowitą energię układu. Zakładając, że siły tarcia są zaniedbywane, całkowita energia jest zachowywana przez system analityczny. Liczbowo jest to zupełnie inna kwestia.

Przed wypróbowaniem podwójnego wahadła wypróbuj prosty wahadło. Zauważysz, że dla metod małego rzędu Runge-Kutty energia systemu będzie rosła w symulacji numerycznych, zamiast pozostać stała (tak dzieje się w twoich symulacjach: ruch jest niczym z niczego). Aby temu zapobiec, można zastosować metody RK wyższego rzędu (ode45 jest rzędu 4; RK rzędu 8 działałby lepiej). Istnieją inne metody zwane „metodami symplektycznymi”, które są zaprojektowane tak, aby symulacje numeryczne oszczędzały energię. Zasadniczo należy przerwać symulację, gdy tylko energia znacznie wzrośnie w porównaniu do inicjalizacji.

I spróbuj zrozumieć proste wahadło, zanim przejdziesz do podwójnego.