Numeryczne rozwiązanie równań różniczkowych

W procesie tworzenia symulacji fizycznych zaczynamy od wymyślenia modelu matematycznego i znalezienia równań różniczkowych, które opisują zależności fizyczne. Następnym krokiem jest nakłonienie komputera do rozwiązania równań, w procesie o nazwie analiza numeryczna.

Rozwiązanie analityczne

W przypadku prostych modeli można użyć rachunku różniczkowego, trygonometrii i innych technik matematycznych, aby znaleźć funkcję, która jest dokładnym rozwiązaniem równania różniczkowego. Nazywa się to rozwiązaniem analitycznym (ponieważ używasz analizy, aby to rozgryźć). Jest to również określane jako rozwiązanie w postaci jawnej.

Preferowane jest rozwiązanie analityczne, ponieważ można uzyskać pewien wgląd w przebieg procesu na podstawie równania. Na przykład w analitycznym rozwiązaniu symulacji pojedynczej sprężyny możemy zobaczyć, jak masa m lub współczynnik sprężystości sprężyny k wpływa na częstotliwość drgań, po prostu patrząc na równanie:

$$x(t) = x_0 \cos(\sqrt{k/m} \; t)$$

W przeciwieństwie do tego, przy rozwiązaniu numerycznym trzeba by przeprowadzić wiele eksperymentów, aby ustalić takie związki; tak jak w prawdziwym życiu, trzeba by eksperymentować z różnymi masami lub sprężynami, dokonywać pomiarów, a następnie tworzyć wykres pokazujący, jak częstotliwość zmienia się dla różnych mas lub wspólczynników sprężystości sprężyn.

Rozwiązanie numeryczne

Większość symulacji fizycznych jest zbyt skomplikowana, aby znaleźć rozwiązanie analityczne. Zamiast tego stosujemy analizę numeryczną, żeby znaleźć przybliżone rozwiązanie numeryczne, które jest tylko szeregiem liczb. Te liczby są wartościami zmiennych modelu w kolejnych momentach czasu. Oto na przykład kilka pierwszych momentów symulacji pojedynczej sprężyny

      t          x       v
  0.00000  −2.00000  0.00000
  0.02500  −1.99626  0.29906
  0.05000  −1.98507  0.59552
  0.07500  −1.96651  0.88827
  0.10000  −1.94070  1.17623
  0.12500  −1.90775  1.45837
  0.15000  −1.86783  1.73364
  0.17500  −1.82113  2.00105
  0.20000  −1.76785  2.25965
  0.22500  −1.70823  2.50851
  0.25000  −1.64252  2.74675

Ta lista liczb jest tym, co napędza wizualną reprezentację symulacji: dla bieżącej wartości czasu t przesuwamy blok do położenia danego przez zmienną x po tym czasie.

Jak to działa? Jak wygenerować tę listę liczb?

Zaczynamy od początkowych wartości zmiennych, które są w pierwszym wierszu powyższej tabeli. Następnie używamy równań różniczkowych dla tej symulacji do obliczania wartości zmiennych po bardzo krótkim czasie. Równanie różniczkowe mówi nam o szybkości zmian każdej zmiennej w dowolnym momencie. Dodając do siebie małe zmiany w czasie (znane również jako całkowanie po czasie) możemy przewidzieć przyszłość.

Równania różniczkowe dla symulacji pojedynczej sprężyny to:

   x' = v (1)
   v' = −km xbm v (2)

Wartości parametrów, takich jak współczynnik sprężystości k , masa m , stała tłumienia b wynoszą odpowiednio:

k = 3.0
m = 0.5
b = 0.1

Przy tych parametrach równania (1) i (2) przyjmują postać

   x' = v (3)
   v' = −6 x − 0.2 v (4)

Warunki początkowe w momencie t = 0 to

x = −2.0
v = 0.0

Możemy znaleźć szybkość zmian zmiennych w momencie t = 0 podstawiając powyższe wartości do równań (3) i (4)

x' = 0
v' = −6 (−2) − 0.2 (0) = +12

Ponieważ znamy początkowe wartości zmiennych i ich szybkość zmian, możemy przewidzieć przyszłe wartości. Przynajmniej przez bardzo krótki okres czasu w przyszłości, przy założeniu, że ta szybkość nie zmienia się zbytnio!

Metoda Eulera obliczania równań różniczkowych

Aby pokazać, jak działa całkowanie numeryczne, użyjemy najprostszej możliwej metody: metody Eulera. Jest prosta, ale niezbyt dokładna. Nieco poniżej wymienimy dokładniejsze metody.

Przy rozwiązywaniu równań różniczkowych (3) i (4), załóżmy, że znamy stan układu po czasie tn określony przez położenie xn i prędkość vn .

Aby uzyskać stan układu w momencie tn+1 = tn + Δt po upływie małego czasu Δt :

   xn+1 = xn + Δt vn (5)
   vn+1 = vn + Δt (− 6 xn − 0.2 vn) (6)

Równania (5) i (6) wykorzystują szybkości zmian dane przez równania różniczkowe (3) i (4), aby uzyskać przybliżone oszacowanie nowej wartości po upływie małego czasu Δt .

Teraz możemy użyć równań (5) i (6), aby dokonać przybliżonych prognoz. Zaczynamy w momencie t = 0 z

x0 = −2
v0 = 0

Przy n = 0 i małym kroku czasowym Δt = 0.025 , równania (5) i (6) dają stan po czasie t1 = t0t = 0.025

x1 = x0 + 0.025 v0 = −2 + 0.025 * 0 = −2
v1 = v0 + 0.025 (−6 x0 − 0.2 v0) = 0 + 0.025 (−6 (−2) − 0.2 (0)) = 0.3

Przy n = 1 w równaniach (5) i (6) możemy uzyskać stan po czasie t2 = t1 + Δt = 0.050

x2 = x1 + 0.025 v1 = −2 + 0.025 (0.3) = −1.9925
v2 = v1 + 0.025 (−6 x1 − 0.2 v1) = 0.3 + 0.025 (−6 (−2) − 0.2 (0.3)) = 0.5985

Przy n = 2 w równaniach (5) i (6) możemy uzyskać stan po czasie t3 = t2 + Δt = 0.075

x3 = x2 + 0.025 v2 = −1.9925 + 0.025 (0.5985) = −1.9775375
v3 = v2 + 0.025 (−6 x2 − 0.2 v2) = 0.5985 + 0.025 (−6 (−1.9925) − 0.2 (0.5985)) = 0.8943825

Nasze przybliżone rozwiązanie numeryczne do tej pory:

      t         x        v
  0.00000  −2.00000  0.00000
  0.02500  −2.00000  0.30000
  0.05000  −1.99250  0.59850
  0.07500  −1.97754  0.89438

Tutaj ponownie dla porównania mamy rozwiązanie, które jest obliczane w symulacji pojedyncza sprężyna, przy użyciu bardzo dokładnej metody Rungego Kutty

      t         x        v
  0.00000  −2.00000  0.00000
  0.02500  −1.99626  0.29906
  0.05000  −1.98507  0.59552
  0.07500  −1.96651  0.88827

Nasze uproszczone rozwiązanie wydaje się niezbyt odległe! Ale jeśli porównamy to z rzeczywistą masą na sprężynie, zobaczymy, że ta metoda jest bardzo niedokładna. Także jeśli spojrzymy na to, czy zachowana jest energia, również zobaczymy niską dokładność tej metody. W przypadku wielu symulacji metoda ta jest tak zła, że staje się niestabilna i „wybucha” z rozwiązaniami dążącymi do nieskończoności.

Inne metody numeryczne

Istnieje kilka metod numerycznego rozwiązywania równań różniczkowych. Mają różne poziomy kompromisu dokładności z kosztami obliczeniowymi. Artykuł w Wikipedii (en) Numerical methods for ordinary differential equations opisuje kilka z nich.

Metody wykorzystywane w symulacjach myphysicslab to:

Te metody numeryczne są definiowane za pomocą równań, które mogą być nieco dezorientujące. Aby przyzwyczaić się do notacji, oto jak można opisać prostą metodę Eulera, której użyliśmy powyżej. W bardziej ogólnej formie możemy zapisać równania (5) i (6) w ten sposób:

   xn+1 = xn + Δt  f (xn, vn) (7)
   vn+1 = vn + Δt  g(xn, vn) (8)

gdzie funkcje f, g z rozwiązywanych przez nas równań różniczkowych są zapisywane tak:

   x' = f (x, v) (9)
   v' = g(x, v) (10)

Problem z prostą metodą Eulera polega na tym, że wykorzystuje ona tylko szybkość zmian na początku przedziału czasu od tn do tn + Δt . A przecież ta szybkość zmian (pochodna) sama się zmienia w tym okresie! Dokładniejsze metody numeryczne szacują, w jaki sposób pochodne w równaniach (9) i (10) zmieniają się w przedziale czasowym.

Źródło: