Na tych laboratoriach skupimy się na scałkowaniu równania ruchu: \[ \mathbf{M} \ddot{\mathbf{x}} = \mathbf{f} - \mathbf{S} \mathbf{x} \] Gdzie \(\mathbf{x}\) to odkształcenie, \(\mathbf{f}\) to siła, \(\mathbf{M}\) to macierz masowa, zaś \(\mathbf{S}\) to macierz sztywności.
Na początek przez \(\mathbf{y}\) oznaczmy prędkość odkształcenia, czyli \(\mathbf{y} = \dot{\mathbf{x}}\). Teraz mamy układ równań pierwszego rzędu:
\[ \begin{cases} \mathbf{M} \dot{\mathbf{y}} & = \mathbf{f} - \mathbf{S} \mathbf{x} \\ \dot{\mathbf{x}} & = \mathbf{y} \\ \end{cases} \]
Zastępując pochodną po lewej stronie przez różnicę skończoną mamy:
\[ \begin{cases} \mathbf{M} \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\text{dt}} & = \mathbf{f} - \mathbf{S} \mathbf{x} \\ \frac{\mathbf{x}^{n+1} - \mathbf{x}^n}{\text{dt}} & = \mathbf{y} \\ \end{cases} \]
Po prawej stronie równania możemy użyć \(\mathbf{x}\) i \(\mathbf{y}\) z nowej (\(n+1\)), bądź starej (\(n\)) iteracji. W zależności co użyjemy otrzymamy mniej lub bardziej uwikłane równanie, a schemat będzie jawny (explicit) bądź niejawny (implicit).
Uwaga: aby porównać różne schematy, każdy schemat zapisz w oddzielnej funkcji, której nagłówek powinien mieć postać:
void Dynamics_schemat(int N, double *x, double *y, double *f, double t_max, double dt);gdzie x i y to początkowe wartości \(\mathbf{x}\) i \(\mathbf{y}\), f to wektor sił, t_max to całkowity czas całkowania, a dt to krok czasowy.
Na początek wstawmy po prawej stronie wartości ze starej iteracji. Otrzymamy: \[ \begin{cases} \mathbf{M} \mathbf{y}^{n+1} & = \mathbf{M} \mathbf{y}^n + \text{dt}(\mathbf{f} - \mathbf{S} \mathbf{x}^n) \\ \mathbf{x}^{n+1} & = \mathbf{x}^n + \text{dt} \mathbf{y}^n \\ \end{cases} \]
Napisz funkcję mnożącą przez macierz masową \(\mathbf{M}\). W pliku MesLib.h jest ona zdefiniowana w analogiczny sposób jak macierz sztywności: przez globalną tablicę M i globalną stałą Mm.
Uwaga: W mnożeniu przez macierz masową, należy także zamrozić wybrane stopnie swobody.
Napisz funkcję całkującą równanie ruchu układu według następującego schematu:
Przeanalizuj dla jakich \(\text{dt}\) układ jest stabilny, a dla jakich nie.
Jak wygląda wzór na całkowitą energię układu (energia potencjalna sprężystości + energia kinetyczna)? Zróżniczkuj ją po \(t\) i pokaż, że jest stała.
Wydrukuj w konsoli jak zmienia się całkowita energia układu w czasie.
Prostą modyfikacją jest użycie po prawej stronie \(\mathbf{x}\) ze starej iteracji i \(\mathbf{y}\) z nowej, otrzymując: \[ \begin{cases} \mathbf{M} \mathbf{y}^{n+1} & = \mathbf{M} \mathbf{y}^n + \text{dt}(\mathbf{f} - \mathbf{S} \mathbf{x}^n) \\ \mathbf{x}^{n+1} & = \mathbf{x}^n + \text{dt} \mathbf{y}^{n+1} \\ \end{cases} \]
Przekopiuj funkcję Dynamics i zmodyfikuj ją tak aby układ na \(\mathbf{y}^{n+1}\) był rozwiązywany przed obliczeniem \(\mathbf{x}^{n+1}\).
Przeanalizuj dla jakich \(\text{dt}\) układ jest stabilny. Wydrukuj zmienność energii.
Możemy także po prawej stronie wziąć obie wartości z nowej iteracji, otrzymując: \[ \begin{cases} \mathbf{M} \mathbf{y}^{n+1} & = \mathbf{M} \mathbf{y}^n + \text{dt}(\mathbf{f} - \mathbf{S} \mathbf{x}^{n+1}) \\ \mathbf{x}^{n+1} & = \mathbf{x}^n + \text{dt} \mathbf{y}^{n+1} \\ \end{cases} \]
Wstawiając drugie równanie do pierwszego otrzymujemy: \[ \mathbf{M} \mathbf{y}^{n+1} = \mathbf{M} \mathbf{y}^n + \text{dt}(\mathbf{f} - \mathbf{S} (\mathbf{x}^n + \text{dt} \mathbf{y}^{n+1})) \] Po przekształceniu: \[ (\mathbf{M} + \text{dt}^2\mathbf{S}) \mathbf{y}^{n+1} = \mathbf{M} \mathbf{y}^n + \text{dt}(\mathbf{f} - \mathbf{S} \mathbf{x}^n) \]
Napisz funkcję mnożącą przez \(\mathbf{M} + \text{dt}^2 \mathbf{S}\)
Zmodyfikuj kod, by realizował schemat w pełni niejawny, zamieniając macierz \(\mathbf{M}\) na \(\mathbf{M} + \text{dt}^2 \mathbf{S}\) w obliczeniu \(y\)-ka.
Przeanalizuj dla jakich \(\text{dt}\) układ jest stabilny. Wydrukuj zmienność energii.
Ostatnia z metod, którymi się zajmiemy bierze po prawej stronie średnią z wartości w nowej i starej iteracji: \[ \begin{cases} \mathbf{M} \mathbf{y}^{n+1} & = \mathbf{M} \mathbf{y}^n + \text{dt} \left( \mathbf{f} - \mathbf{S} \frac{\mathbf{x}^{n+1} + \mathbf{x}^n}{2} \right) \\ \mathbf{x}^{n+1} & = \mathbf{x}^n + \text{dt} \frac{\mathbf{y}^{n+1} + \mathbf{y}^n}{2} \\ \end{cases} \]
Po wstawieniu drugiego równania do pierwszego mamy: \[ \mathbf{M} \mathbf{y}^{n+1} = \mathbf{M} \mathbf{y}^n + \text{dt} \left( \mathbf{f} - \mathbf{S} \frac{\mathbf{x}^n + \text{dt} \frac{\mathbf{y}^{n+1} + \mathbf{y}^n}{2} + \mathbf{x}^n}{2} \right) \] Po uproszczeniu: \[ \mathbf{M} \mathbf{y}^{n+1} = \mathbf{M} \mathbf{y}^n + \text{dt} \left( \mathbf{f} - \mathbf{S} \left( \mathbf{x}^n + \text{dt} \frac{\mathbf{y}^{n+1} + \mathbf{y}^n}{4} \right) \right) \] Ostatecznie: \[ \left( \mathbf{M} + \frac{\text{dt}^2}{4} \mathbf{S} \right) \mathbf{y}^{n+1} = \mathbf{M} \mathbf{y}^n + \text{dt} \left( \mathbf{f} - \mathbf{S} \left( \mathbf{x}^n + \text{dt} \frac{\mathbf{y}^n}{4} \right) \right) \]
Napisz funkcję mnożącą przez \(\mathbf{M} + \frac{\text{dt}^2}{4} \mathbf{S}\).
Napisz funkcję całkującą równanie ruchu według następującego schematu (dla uproszczenia zapisu pominęliśmy number n iteracji):
Przeanalizuj dla jakich \(\text{dt}\) układ jest stabilny. Wydrukuj zmienność energii.
Udowodnij, że metoda ,,w pół kroku’’ zachowuje energię całkowitą układu.1
Podpowiedź: tak jak \(a^2-b^2=(a+b)(a-b)\) to \(x_{n+1}^TMx_{n+1} - x_{n}^TMx_{n} = (x_{n+1} - x_{n})^TM(x_{n+1} + x_{n})\)↩︎