Wdrażanie pół-niejawnego Eulera wstecznego w systemie 1-DOF z masą sprężynową

Question

Mam prosty (masowy) system z dwoma punktami połączonymi ze sprężyną. Jeden punkt jest ustawiony na suficie, więc chcę obliczyć położenie drugiego punktu za pomocą metody numerycznej. Tak więc, zasadniczo dostaję pozycję drugiego punktu i jego prędkość, i chcę wiedzieć, jak te dwie wartości aktualizują się po jednym taktowaniu.

następujące siły ma wpływ na punkcie:

• siły grawitacji, przeprowadzonego przez -g * m
• siły sprężyny, podanych k * (l - L) k jest sztywność, L jest obecny długość i L oznacza początkowy długość
• siłę tłumiącą, ponieważ przez -d * v

zsumowane, co prowadzi do

• F = -g * m + k * (l - L)
• Fd = -d * v

Stosując na przykład Explicit Euler, można wyprowadzić następujące:

• newVelocity = oldVelocity + dt * (F + Fd)/m, używając F = m * a.

Jednak teraz chcę używać wpół ukryte tyłu Euler, ale nie można dokładnie dowiedzieć się, gdzie wyprowadzić Jacobians z itp

Answer 1

Więc to chyba najłatwiej zobaczyć, jak to idzie z najpierw należy rozważyć metodę całkowicie niejawną, a następnie przejść do niejako pośrednio.

niejawny Euler musiałaby (nazwijmy te eqn (1)):

newPos = oldPos + dt * newVelocity 
newVelocity = oldVelocity + dt * (-g * m + k*(newPos - L) - d*newVelocity)/m 

Na razie niech po prostu zmierzyć pozycje w stosunku do L, dzięki czemu możemy pozbyć się tego -KL perspektywie. Przekształcając możemy skończyć z

(newPos, newVelocity) - dt * (newVelocity, k/m newPos - d/m newVelocity) = (oldPos, oldVelocity - g*dt) 

i oddanie, które w postaci macierzowej

((1,-dt),(k/m, 1 - d/m)).(newPos, newVelocity) = (oldPos, oldVelocity -g*dt) 

Jeżeli wiesz wszystko na matrycy, a wszystko na RHS, i po prostu trzeba rozwiązać dla wektora (newPos, newVelocity). Możesz to zrobić przy pomocy dowolnego rozwiązania Ax = b (eliminacja gaussowska ręcznie działa w tym prostym przypadku). Ale ponieważ wspomniałeś o Jacobianach, prawdopodobnie chcesz rozwiązać ten problem, używając iteracji Newtona-Raphsona lub czegoś podobnego.

W takim przypadku, jesteś w istocie chce rozwiązać zera równania

((1,-dt),(k/m, 1-d/m)).(newPos, newVelocity) - (oldPos, oldVelocity -g*dt) = 0 

co znaczy, F (newPos, newVelocity) = (0,0).Masz poprzednią wartość do użycia jako początkowe odgadnięcie, (oldPos, oldVelocity). Teraz po prostu chcesz iteracyjne na

(x, v) ^{n + 1} = (x, v) ⁿ + f ((x, v) ⁿ)/f '((x, v) ⁿ)

, dopóki nie uzyskasz wystarczająco dobrej odpowiedzi. Tutaj

f(newPos,newVel) = ((1,-dt),(k/m, 1-d/m)).(newPos, newVelocity) - (oldPos, oldVelocity -g*dt) 

oraz f "(newPos, newVel) jest jakobian odpowiada matrycy

((1,-dt),(k/m, 1-d/m)) 

przechodzi proces dla pół-niejawny jest taka sama, ale trochę łatwiej - nie wszyscy terminy RHS w eqns (1) to nowe ilości. Jak to się zwykle robi to

newPos = oldPos + dt * newVelocity 
newVelocity = oldVelocity + dt * (-g * m + k*oldPos - d*newVelocity)/m 

np prędkość zależy od starej wartości czasowej pozycji i pozycji na nowej wartości czasowej prędkości. (Jest to bardzo podobne do integracji "leapfrog"). Powinieneś być w stanie wykonać powyższe kroki całkiem łatwo dzięki temu nieco odmiennemu zestawowi równań. Zasadniczo, k/m termin w macierzy powyżej spada.