forked from ccfd/courses
-
Notifications
You must be signed in to change notification settings - Fork 0
/
info2_lab05.Rmd
157 lines (131 loc) · 9.12 KB
/
info2_lab05.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
---
number: 5
course: Informatyka 2
material: Instrukcja 5
author: W. Gryglas
---
# Całkowanie numeryczne układów równań
Pliki do wykorzystania w poniższym ćwiczenie można pobrać za pomocą poniższych linków:
- [Plik nagłówkowy rk4.h](http://ccfd.github.io/courses/code/info2/rk4.h)
- [Plik źródłowy rk4.cpp](http://ccfd.github.io/courses/code/info2/rk4.cpp)
## 1. Wstęp
Z poprzednich ćwiczeń dowiedzieliśmy się jak można dokonać całkowania równań różniczkowych zwyczajnych (z ang. Ordinary Differential Equations, ODE). Opisane metody pozwalają na znalezienie konkretnej funkcji spełniającej dowolne równanie pierwszego rzędu, tj. postaci:
$$ \frac{dy}{dt} = f(t, y(t)) $$
Jednak świat który opisujemy za pomocą równań różniczkowych nie jest niestety tak prosty. W wielu zagadnieniach napotykamy na sytuacje gdy należy "jednocześnie" rozwiązać dwa równania różniczkowe. Mówiąc "jednocześnie" mamy tutaj na myśli sytuację, w której w obu równaniach pojawiają się nieznane(których poszukujemy rozwiązując równania różniczkowe) funkcje, np.:
$$
\begin{aligned}
\frac{dx}{dt} &= sin(y) + cos(x) \cdot t \\
\frac{dy}{dt} &= sin(y) \cdot t + cos(x)
\end{aligned}
$$
Jak widać w powyższym układzie równań poszukiwane funkcje to $x = x(t)$ i $y=y(t)$. Obie są funkcjami zmiennej niezależnej $t$. W obu równaniach pojawia się jednak zarówno funkcja $x$ jak i $y$. Równania te są sprzężone, a zatem nie da się w ogólności rozwiązać (scałkować) tego układu w sposób skwencyjny, tj. znaleźć najpierw pierwszą funkcję, a potem oddzielnie drugą. W szczególnych przypadkach można się pokusić o transofrmację wyjściowych równań, tak aby wyrugować jedną ze zmiennych. Jednak w ogólnym przypadku, gdy mamy doczynienia z silnie uwikłanymi funkcjami prawych stron, taki zabieg jest nie możliwy. Dlatego aby znaleźć poszukiwane rozwiązania należy całkować równocześnie wszystkie równania. Takim zagadnieniem zajemiemy się w poniższym ćwiczeniu. Postaramy się zastosować poznane metody z poprzednich zajęć do rozwiązania układów równań pierwszego rzędu.
Problem rozwiązywania układów rówań pierwszego rzędu bardzo często pojawia się w związku z problemami opisanymi równaniami różniczkowymi zwyczajnymi wyżych rzędów. Równanie różniczkowe zwyczajne drugiego rzędu są jednymi z bardziej rozpowszechnionych równań w mechanice, ponieważ służą do opisu dynamiki, np. ruchu bryły sztywnej pod działaniem sił zewnętrznych. Równanie te przyjmują następującą postać:
$$ \frac{d^2x}{dt^2} = f(t, x, \frac{dx}{dt}) $$
Każde tego typu równanie można przedstawić w równoważnej formie układu rówań **pierwszego** rzędu:
$$
\begin{aligned}
\frac{dv_x}{dt} &= f_1(t, x, v_x) \\
\frac{dx}{dt} &= v_x
\end{aligned}
$$
Powyższy układ jest tożsamy z równanie drugiego rzędu i znalezienie przebiegu funkcji $x=x(t)$ będzie odpowiadało scałkowaniu rówania drugiego rzędu. Tego typu problemem zajmiemy się w drugiej części tego ćwiczenia.
## 2. Całkowanie układów rówań
Na poprzednich zajęciach poznaliśmy dwie metody jawnego całkowania równań różniczkowych zwyczajnych pierwszego rzędu:
- Eulera
- Rungego-Kutty
Metody te należą do klasy metod iteracyjnych i pozwalają znaleź przebieg poszukiwanej funkcji obliczając jej wartość w kolejnych krokach czasowych. W przypadku układów równań czynimy podbnie, tj. obliczamy nowe wartości funkcji bazując na poprzednim rozwiązaniu. Różnica w tym przypadku polega na tym, że musimy kolejny krok całkowania wykonać dla wszystkich równań układu. Weźmy pod uwagę przykładowy układ równań ze wstępu do tych ćwiczeń:
$$
\begin{aligned}
\frac{dx}{dt} &= sin(y) + cos(x) \cdot t \\
\frac{dy}{dt} &= sin(y) \cdot t + cos(x)
\end{aligned}
$$
z warunkami początkowymi:
$$
x(t=0) = 0 \\
y(t=0) = 0
$$
Dla tagiego układu równań całkowanie metodą Eulera można opisać za pomocą następującego algorytmu:
1. Ustaw $x_0 = x(t=0)$, $y_0 = y(t=0)$, wartość czasu $t=0$ oraz wybierz krok całkowania $dt$, np. $dt = 0.001$.
2. Oblicz pochodne funkcji $x$ i $y$ wykorzystując znane wielkości $x_0, y_0, t$, w tym przypadku:
$$
dxdt = sin(y_0) + cos(x_0) \cdot t = sin(0) + cos(0) \cdot 0 = 0 \\
dydt = sin(y_0) \cdot 0 + cos(x_0) = sin(0) \cdot 0 + cos(0) = 1
$$
3. Wykorzystując wartości pochodnych dla czasu $t=0$ oblicz nowe wartości funkcji $x$, $y$:
$$
x = x_0 + dt \cdot dxdt = 0 + 0.001 \cdot 0 = 0 \\
y = y_0 + dt \cdot dydt = 0 + 0.001 \cdot 1 = 0.001
$$
4. Zaktualizuj zmienne podstawiając:
$$
x_0 = x = 0
y_0 = y = 0.001
t = t + dh = 0.001
$$
5. Przejdź do punktu 2 i potórz operacje aż zmienna $t$ osiągnie wartość końcową $t_k$
Korzystając z powyższego algorytmu można uzyskać następujący wynik dla czasu w zakresie [0,1]
```{r, echo=FALSE}
x_0 = 0; y_0 = 0; dt = 0.001; t = 0
X=x_0; Y=y_0
for(i in 1:1000)
{
X = append(X, x_0 + dt*(sin(y_0) + cos(x_0)*t) )
Y = append(Y, y_0 + dt*(sin(y_0)*t + cos(x_0)) )
x_0 = X[i]; y_0 = Y[i]
t = t + dt
}
t = seq(from=0, to=1, by=dt)
plot(t, X, type="l", xlab="t",ylab="Wartosci", col="blue")
lines(t, Y, col="red")
legend(0.05,y=0.3, legend = c("x(t)","y(t"), col=c("blue","red"), lty=c(1,1), cex=1.2)
```
Powyższy algorytm bedzie bardzo podobny dla metody Rungego-Kutty. Różnica będzie jedynie polegała na zastąpieniu sposobu obliczania wartości funkcji parawych stron, tj. w tym przypadku zmiennych $dxdt$ oraz $dydt$.
### Ćwiczenia
1. Zaimplementuj opisany algorytm i dokonaj całkowania układu równań przedstawionego powyżej.
2. Wartości rozwiązania dla kolejnych kroków zapisz w tablicach X, Y.
3. Wyświetl wartości funkcji X(t) oraz Y(t) na spólnym wykresie używając biblioteki graficznej.
4. Zmodyfikuj program tak aby całkowanie odbywało się przy wykorzystaniu metody Rungego-Kutty.
## 3. Sprowadzanie równań wyższego rzędu do układów równań
Teraz przejdźmy do kolejnego zagadnienia. Aby zastosować poznane metody do problemów opisanych równaniami wyższych rzędów należy przekształcić wyjściowy problem do układu równań pierwszego rzędu. Aby przeanalizować ten proces weźmy pod uwagę następujące zagadnienie:
Ruch wózka, zakładając, że w pozycji wyjściowej sprężyna jest nierozciągnięta oraz masa kół jest pomijalna, można opisać za pomocą prostego równania dynamiki bazując na III zasadzie dynamiki Newtona:
$$ F_{wózek} = - k\cdot x = m \cdot a$$
wiedząc, że przyśpieszenie to druga pochodna przemieszczenia po czasie możemy ostatecznie zapisać:
$$ \frac{d^2x}{dt^2} = - \frac{k}{m} \cdot x $$
Oprócz równania potrzebujemy jeszcze warunków początkowych. Z racji tego, że jest to równanie drugiego rzędu będziemy potrzebowali dwóch takich warunków. W tym przypadku możemy wprost założyć, że startowe położenie wózka jest w początku układu współrzędnych, czyli:
$$ x(t=0) = 0 $$
oraz, że na początku wózek jest w spoczynku:
$$ \frac{dx}{dt}(t=0)=0 $$
Teraz musimy przekształcić ten problem do problemu opisanego układem dwóch równań różniczkowych pierwszego rzędu. Można to uczynić na kilka sposobów. Jednym z nich jest skorzystanie z równań kanonicznych, które opsiują ruch ciała za pomocą pary [pęd, położenie]. Innym podjeściem, które najłatwiej w tej chwili zrealizować, bedzie przekształcenie równania wyjściowego przez wprowadzenie nowej zmiennej - $v$ i zastosowanie prostego podstawienia:
$$ v = \frac{dx}{dt} $$
Powyższe podstawienie automatycznie definiuje nam nowe równanie, które wiąże prędkość i położenie. Teraz, pozostaje tylko zastosować powyższe podstawienie w kontekście równania wyjściowego:
$$ \frac{d^2x}{dt^2} = \frac{d}{dt} (\frac{dx}{dt}) = \frac{dv}{dt} = -\frac{k}{m}x $$
Oprócz równań podobnie należy uczynić z warunkami początkowymi. Ostatcznie otrzymujemy następujący układ równań:
$$
\frac{dv}{dt} = -\frac{k}{m} x \\
\frac{dx}{dt} = v \\
x(t=0) = 0 \\
v(t=0) = 0
$$
### Ćwiczenia
Ruch wahadła matematycznego opisuje równanie różniczkowe z warunkami początkowymi (pomijamy
opór powietrza, nić jest nieważka i nierozciągliwa):
$$
\begin{aligned}
\frac{d^2\alpha}{dt^2} &= -\frac{g}{l}sin(\alpha) \\
\alpha(t_0) &= \alpha_0 \\
\frac{d\alpha}{dt}(t_0) &= \omega_0
\end{aligned}
$$
gdzie:
$\alpha$ - kąt wychylenia wahadła z położenia równowagi,
$g$ - przyśpieszenie ziemskie,
$l$ - długość wahadła,
$m$ - masa kuli zaczepionej na końcu wahadła.
1. Sprowadzić powyższe równanie różniczkowe rzędu drugiego do układu równań różniczkowych rzędu pierwszego.
2. Napisać program który używając metody Rungego-Kutty wyznacza zależność kąta wychylenia wahadła od czasu $\alpha(t)$ oraz prędkość kątową $\omega = \frac{d\alpha}{dt}$ dla czasu $t \in \left[0,10\right]$
3. Dodatkowo wyznaczyć zależność energii całkowitej wahadła od czasu $E(t)$, energia całkowita wahadła wyraża się wzorem:
$$E = \frac{ml^2}{2}(\frac{d\alpha}{dt})^2 + mgl(1-cos(\alpha))$$
**Uwaga**: Przy braku dyssypacji, energia mechaniczna powinna być stała.
4. Powtórzyć obliczenia dla różnych kroków czasowych.