Autor: A
. A
aniewski-Wołłk Creative Commons License:
MetNum: Lab 5
Attribution Share Alike
Metody Numeryczne: Instrukcja 5 Zadanie 2 Napisz funkcję całkującą równanie ruchu układu wg. następującego
schematu:
Całkowanie równania ruchu
" Oblicz b = My + dt(F - Sx)
Na tych laboratoriach skupimy się na scałkowaniu równania ruchu:
" Oblicz x = x + dty
Mć = F - Sx
" Rozwiąż układ: My = b
Gdzie x to odkształcenie, M to macierz masowa, zaś S to macierz sztywności.
" Co 10-tą iterację wyświetl belkę.
Na poczÄ…tek przez y oznaczmy prÄ™dkość odksztaÅ‚cenia, czyli y = ‹. Teraz
Zadanie 3 Przeanalizuj dla jakich dt układ jest stabilny, a dla jakich nie.
mamy układ równań pierwszego rzędu:
Zadanie 4 Jak wygląda wzór na całkowitą energię układu (energia potencjalna
MŹ = F - Sx
sprężystości + praca sił + energia kinetyczna)? Zróżniczkuj ją po t i pokaż, że
‹ = y
jest stała
Zastępując pochodną po lewej stronie przez różnice skończoną mamy:
Zadanie 5 Wydrukuj w konsoli jak zmienia się całkowita energia układu w
czasie.
yn+1-yn
M = F - Sx
dt
xn+1-xn
= y
dt
2 Schemat pół niejawny (semi-implicit)
Po prawej stronie równania możemy użyć x i y z nowej (n + 1), bądz
starej
(n) iteracji. W zależności co użyjemy otrzymamy mniej lub bardziej uwikłane
Prostą modyfikacją jest użycie po prawej stronie x ze starej iteracji i y z nowej,
równanie, a schemat będzie jawny (explicit) bądz
niejawny (implicit).
otrzymujÄ…c:
Uwaga: by porównać różne schematy, każdy schemat napisz w nowej funkcji,
Myn+1 = Myn + dt(F - Sxn)
która powinna brać następujące argumenty:void Dynamics(int n, double
xn+1 = xn + dtyn+1
* x, double * y, double T, double dt);, gdziexiyto poczÄ…tkowe war-
tości x i y,Tto całkowity czas całkowania, adtto krok czasowy.
Zadanie 6 Zmodyfikuj kod rozwiązując układ na y przed modyfikacją x-a.
Zadanie 7 Przeanalizuj dla jakich dt układ jest stabilny. Wydrukuj zmienność
1 Schemat prawie jawny (almost explicit)
energii.
Na początek wstawmy po prawej stronie wartości ze starej iteracji. Otrzymamy:
3 Schemat niejawny (fully-implicit)
Myn+1 = Myn + dt(F - Sxn)
Możemy także po prawej stronie wziąć obie wartości z nowej iteracji, otrzymu-
xn+1 = xn + dtyn
jÄ…c:
Zadanie 1 Napisz funkcję mnożącą przez macierz masową. W plikuMesLib.h
Myn+1 = Myn + dt(F - Sxn+1)
jest ona zdefiniowana w analogiczny sposób jak macierz sztywności: przez ma-
xn+1 = xn + dtyn+1
cierz M i stałą Mm. UWAGA: W mnożeniu przez macierz masową, należy
Wstawiając drugie równanie do pierwszego otrzymujemy:
także zamrozić wybrane stopnie swobody.
Myn+1 = Myn + dt(F - S(xn + dtyn+1))
Wydział Mechaniczny Energetyki i Lotnictwa, Politechnika Warszawska 1
Autor: A
. A
aniewski-Wołłk Creative Commons License:
MetNum: Lab 5
Attribution Share Alike
dt
Przekształcając: " Oblicz x = x + y
4
dt2
(M + dt2S)yn+1 = Myn + dt(F - Sxn)
" Rozwiąż układ: (M + S)y = b
4
dt
Zadanie 8 Napisz funkcję mnożącą przez M + dt2S
" Oblicz x = x + y
2
Zadanie 9 Zmodyfikuj kod, by realizował schemat w pełni niejawny, zamienia- " Co 10-tą iterację wyświetl belkę.
jÄ…c macierz M na M + dt2S w obliczeniu y-ka
Zadanie 13 Przeanalizuj dla jakich dt układ jest stabilny. Wydrukuj zmien-
Zadanie 10 Przeanalizuj dla jakich dt układ jest stabilny. Wydrukuj zmien- ność energii.
ność energii.
Zadanie 14 Udowodnij, że metoda pół kroku zachowuje energię układu.1
4 W pół kroku (midpoint)
Ostatnia z omówionych metod bierze po prawej stronie średnią z wartości w
nowej i starej iteracji:
xn+1+xn
Myn+1 = Myn + dt(F - S )
2
yn+1+yn
xn+1 = xn + dt
2
Po wstawieniu drugiego równania do pierwszego mamy:
yn+1+yn
xn + dt + xn
2
Myn+1 = Myn + dt(F - S )
2
Przekształcając:
yn+1 + yn
Myn+1 = Myn + dt(F - S(xn + dt )
4
Ostatecznie:
dt2 yn
(M + S)yn+1 = Myn + dt(F - S(xn + dt )
4 4
dt2
Zadanie 11 Napisz funkcję mnożącą przez M + S
4
Zadanie 12 Napisz funkcję całkującą równanie ruchu układu wg. następują-
cego schematu:
dt
" Oblicz x = x + y
4
1
Podpowiedz: tak jak a2 - b2 = (a + b)(a - b) to xT Mxn+1 - xT Mxn = (xn+1 -
n+1 n
" Oblicz b = My + dt(F - Sx) xn)T M(xn+1 + xn)
Wydział Mechaniczny Energetyki i Lotnictwa, Politechnika Warszawska 2