Implementazione di Euler semi-implicito a ritroso in un sistema di molle di massa 1-DOF

Question

Ho un semplice sistema (di massa) con due punti collegati a una molla. Un punto è fissato al soffitto, quindi voglio calcolare la posizione del secondo punto usando un metodo numerico. Quindi, in sostanza, ottengo la posizione del secondo punto e la sua velocità, e voglio sapere come questi due valori si aggiornano dopo un timestep.

le seguenti forze hanno effetto sul punto:

• forza gravitazionale, in -g * m
• forza Spring, in k * (l - L) con k essendo la rigidità, l è la lunghezza corrente e L è la iniziale lunghezza
• forza di smorzamento, in -d * v

Riassume, questo porta a

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

applicazione ad esempio di Eulero esplicito, si può ricavare la seguente:

• newVelocity = oldVelocity + dt * (F + Fd)/m, utilizzando F = m * a.

Tuttavia, ora voglio usare semi-implicito Eulero, ma non può esattamente capire dove derivano i Jacobiani da ecc

Answer 1

Quindi è probabilmente più facile vedere come questo va da considerando prima il metodo pienamente implicito, quindi andando al semi-implicito.

implicito di Eulero sarebbe (chiamiamolo questi eqn (1)):

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

Per ora limitiamoci a misurare posizioni relative a L in modo che possiamo sbarazzarci di quel termine -kl. Riordinando si finisce con

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

e la messa che in forma matriciale

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

Dove tu sai tutto nella matrice, e tutto sul RHS, e hai solo bisogno di risolvere per il vettore (newPos, newVelocity). Puoi farlo con qualsiasi Ax = b risolutore (l'eliminazione gaussiana a mano funziona in questo semplice caso). Ma dal momento che menzioni Jacobians, presumibilmente stai cercando di risolvere questo con l'iterazione di Newton-Raphson o qualcosa di simile.

In tal caso, si sta essenzialmente cercando di risolvere gli zeri della equazione

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

vale a dire, f (newPos, newVelocity) = (0,0).Hai un valore precedente da utilizzare come ipotesi iniziale (oldPos, oldVelocity). Ora si vuole solo iterare sulle

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

finché non si ottiene una risposta sufficientemente buona. Qui,

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

e f '(newPos, newVel) è la Jacobiana corrispondente matrice

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

Passando attraverso il processo per semi-implicito è lo stesso, ma un po' più facile - non tutti i termini RHS in eqns (1) sono nuove quantità. Il modo in cui è fatto solitamente è

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

ad esempio, la velocità dipende dal vecchio valore temporale della posizione, e la posizione sul nuovo valore temporale della velocità. (Questo è molto simile all'integrazione "leapfrog".) Dovresti essere in grado di lavorare con i passaggi precedenti piuttosto facilmente con questo insieme leggermente diverso di equazioni. Fondamentalmente, il termine k/m nella matrice sopra si allontana.