matlab — Почему численное решение ОДУ отходит от неустойчивого равновесия?

Почему численное решение ОДУ отходит от неустойчивого равновесия?

Я думаю, что два основных момента уже были высказаны Брайаном и Эртксимом: ваше начальное значение является неустойчивым равновесием, а тот факт, что ваши численные вычисления никогда не бывают точными, обеспечивает небольшое возмущение, которое вызовет нестабильность.

Чтобы немного подробнее описать, как это работает, рассмотрим вашу проблему в форме общей проблемы начальных значений.

begin {уравнение} dot { mathbf {y}} (t) = mathbf {M} ^ {- 1} mathbf {f} (t, mathbf {y} (t)) end {уравнение} где $ mathbf {y} (t) = (x_1 (t), x_2 (t), x_3 (t), x_4 (t)) $ и

begin {уравнение} mathbf {f} (t, mathbf {y} (t)) = begin {bmatrix} x_2 \ -V_1 — G_1 \ x_4 \ -V_2 — G_2 end {bmatrix} конец {} уравнение

В точной арифметике у вас будет $ mathbf {f} (0, mathbf {y} _0) = 0 $ и, следовательно, $ dot { mathbf {y}} (0) = 0 $, и ничего не изменится: ваш Система остается в равновесии. Тем не менее, арифметика в компьютере не является точной, а это означает, что по ряду причин ваша правая часть не точно равна нулю, но равна некоторому $ tilde { mathbf {f}} $, который почти, но не совсем равен нулю.

В своем коде вы можете проверить это, вычисляя

 norm(rhs(0,[pi/2 0 0 0])) 

что дает 6.191e-16 — так почти, но не точно ноль. Как это влияет на динамику вашей системы?

Согласно некоторым предположениям, эффект того, что $ mathbf {f} $ не является точно нулевым, будет таким же, если вы не начнете с начального значения $ mathbf {y} _0 $, которое вы прописываете, но со значения, которое немного отличается давайте назовем это $ tilde { mathbf {y}} _ 0 $ .

Кроме того, в течение очень короткого времени решение вашей системы выглядит как решение линеаризованной системы.

begin {уравнение} dot { mathbf {y}} (t) = mathbf {f} (0, mathbf {y} _0) mathbf {f} ‘(0, mathbf {y} _0) left ( mathbf {y} (t) — mathbf {y} _0 right) = mathbf {f} ‘(0, mathbf {y} _0) left ( mathbf {y} (t) — mathbf {y} _0 right) end {уравнение}

где $ mathbf {f} ‘$ — это якобиан вашей функции $ mathbf {f} $ или rhs в вашем коде. Поскольку $ mathbf {y} _0 $ является константой, мы можем преобразовать это в уравнение для $ mathbf {d} (t): = mathbf {y} (t) — mathbf {y} _0 $ , где $ mathbf {d} $ говорит, как далеко от начального значения мы находимся:

begin {уравнение} dot { mathbf {d}} (t) = mathbf {f} ‘(0, mathbf {y} _0) mathbf {d} (t). конец {} уравнение

Я не мог потрудиться вычислить якобиан вручную, поэтому я использовал автоматическое дифференцирование, чтобы получить хорошее приближение:

begin {уравнение} mathbf {J}: = mathbf {f} ‘(0, mathbf {y} _0) = begin {bmatrix} 0 & 1 & 0 & 0 \ 9.81 & 0 & 2.4525 & 0 \ 0 & 0 & 0 & 1 \ 2.4525 & 0 & 2.4525 & 0 end {bmatrix} end {уравнение}

так что ваше уравнение становится

begin {уравнение} dot { mathbf {d}} (t) = mathbf {J} mathbf {d} (t), mathbf {d} (0) = mathbf { tilde {y}} _0 — mathbf {y} _0 end {уравнение}

Теперь нам нужен один последний шаг: мы можем вычислить разложение по собственным значениям якобиана так, чтобы

begin {уравнение} mathbf {J} = mathbf {Q} mathbf {D} mathbf {Q} ^ {- 1} end {уравнение}

где $ mathbf {D} $ — диагональная матрица, содержащая собственные значения $ mathbf {J} $ и $ mathbf {Q} $ — ортогональные матрицы, представляющие преобразования координат. Затем мы можем преобразовать уравнение для $ mathbf {d} $ в уравнение для $ mathbf {e} (t): = mathbf {Q} ^ {- 1} mathbf {d} (t) $, которое читает

begin {уравнение} dot { mathbf {e}} (t) = mathbf {D} mathbf {e} (t), mathbf {e} (0) = mathbf {Q} ^ {- 1 } mathbf {d} _0. конец {} уравнение

Поскольку $ mathbf {D} $ является диагональю, это фактически четыре независимых уравнения

begin {уравнение} dot {e} _i (t) = lambda_i e_i (t), e_i (0) = i- textrm {th компонент {}} mathbf {Q} ^ {- 1} mathbf { d} _0 end {уравнение}

с $ i = 1,2,3,4 $ . Если вы вычислите собственные значения, вы обнаружите, что наибольшее значение равно $ lambda_1 = 3.2485 $ . Следовательно,

begin {уравнение} e_1 (t) = e_1 (0) e ^ {3.2485 t}. конец {} уравнение

Теперь, если бы арифметика на вашем компьютере была точной, вы бы имели $ mathbf {d} (0) = 0 $ , таким образом, $ mathbf {e} (0) = mathbf {Q} ^ {- 1} mathbf { d} (0) = 0 $ и, следовательно, $ e_1 (0) = 0 $, и ничего не произойдет. Но поскольку это не так, у вас есть небольшой, но конечный $ e_1 (0) $, который экспоненциально усиливается. Отсюда и быстрое отклонение от равновесия в вашем решении.

Понравилась статья? Поделиться с друзьями:
JavaScript & TypeScript
Adblock
detector