Анализ устойчивости фон Неймана

В численный анализ, анализ устойчивости фон Неймана (также известный как Фурье анализ устойчивости) - это процедура, используемая для проверки устойчивости конечно-разностных схем применительно к линейным уравнениям в частных производных. Анализ основан на разложении Фурье числовой ошибки и был разработан в Национальной лаборатории Лос-Аламоса после краткого описания в статье 1947 года автора Британские исследователи Крэнк и Николсон. Этот метод является примером явного интегрирования по времени, где функция, определяющая управляющее уравнение, оценивается в текущее время. Позже этот метод был подвергнут более строгому рассмотрению в статье, написанной в соавторстве с Джоном фон Нейманом.

Числовая устойчивость

Устойчивость числовых схем тесно связана с числовая ошибка. Схема конечных разностей является стабильной, если ошибки, допущенные на одном временном шаге вычисления, не приводят к увеличению ошибок при продолжении вычислений. Нейтрально устойчивая схема - это схема, в которой ошибки остаются постоянными по мере продолжения вычислений. Если ошибки уменьшаются и, в конечном итоге, затухают, численная схема считается стабильной. Если же, наоборот, ошибки растут со временем, численная схема называется неустойчивой. Устойчивость численных схем можно исследовать, выполнив анализ устойчивости по фон Нейману. Для задач, зависящих от времени, стабильность гарантирует, что численный метод дает ограниченное решение, когда решение точного дифференциального уравнения ограничено. Стабильность, в общем, может быть трудно исследовать, особенно когда рассматриваемое уравнение нелинейно.

. В некоторых случаях стабильность по фон Нейману необходима и достаточна для устойчивости в смысле Лакса – Рихтмайера (как используется в Теорема эквивалентности Лакса ): модели PDE и конечно-разностной схемы являются линейными; PDE является постоянным коэффициентом с периодическими граничными условиями и имеет только две независимые переменные; и в схеме используется не более двух временных уровней. Стабильность по фон Нейману необходима в гораздо более широком спектре случаев. Он часто используется вместо более подробного анализа стабильности, чтобы дать хорошее предположение об ограничениях (если таковые имеются) на размеры шага, используемых в схеме, из-за ее относительной простоты.

Иллюстрация метода

Метод фон Неймана основан на разложении ошибок на ряд Фурье. Для иллюстрации процедуры рассмотрим одномерное уравнение теплопроводности

$\partial \; u\; \partial \; t\; =\; \alpha \; \partial \; 2\; u\; \partial \; x\; 2\; \{\backslash \; displaystyle\; \{\backslash \; frac\; \{\backslash \; partial\; u\}\; \{\backslash \; partial\; t\}\}\; =\; \backslash \; alpha\; \{\backslash \; frac\; \{\backslash \; partial\; ^\; \{2\}\; u\}\; \{\backslash \; partial\; x\; ^\; \{2\}\}\}\}$

, определенный на пространственном интервале $L\; \{\backslash \; displaystyle\; L\}$, который может быть дискретизируется как

$(1)\; ujn\; +\; 1\; =\; ujn\; +\; r\; (uj\; +\; 1\; n\; -\; 2\; ujn\; +\; uj\; -\; 1\; n)\; \{\backslash \; displaystyle\; \backslash \; quad\; (1)\; \backslash \; qquad\; u\_\; \{j\}\; ^\; \{n\; +\; 1\}\; =\; u\_\; \{j\}\; ^\; \{n\}\; +\; r\; \backslash \; left\; (u\_\; \{j\; +\; 1\}\; ^\; \{n\}\; -2u\_\; \{j\}\; ^\; \{n\}\; +\; u\_\; \{j-1\}\; ^\; \{n\}\; \backslash \; right)\}$

где

$r\; =\; \alpha \; \Delta \; T\; (\Delta \; x)\; 2\; \{\backslash \; displaystyle\; r\; =\; \{\backslash \; frac\; \{\backslash \; alpha\; \backslash ,\; \backslash \; Delta\; t\}\; \{\backslash \; left\; (\backslash \; Delta\; x\; \backslash \; right)\; ^\; \{2\}\}\}\}$

и решение $ujn\; \{\backslash \; displaystyle\; u\_\; \{j\}\; ^\; \{n\}\}$дискретного уравнения аппроксимирует аналитическое решение $u\; (x,\; t)\; \{\backslash \; displaystyle\; u\; (x,\; t)\}$PDE на сетке.

Определите ошибку округления $ϵ\; jn\; \{\backslash \; displaystyle\; \backslash \; epsilon\; \_\; \{j\}\; ^\; \{n\}\}$как

$ϵ\; jn\; =\; N\; jn\; -\; ujn\; \{\backslash \; displaystyle\; \backslash \; epsilon\; \_\; \{j\}\; ^\; \{n\}\; =\; N\_\; \{j\}\; ^\; \{n\}\; -u\_\; \{j\}\; ^\; \{n\}\}$

где $ujn\; \{\backslash \; displaystyle\; u\_\; \{j\; \}\; ^\; \{n\}\}$- решение дискретизированного уравнения (1), которое было бы вычислено при отсутствии ошибки округления, и $N\; jn\; \{\backslash \; displaystyle\; N\_\; \{j\}\; ^\; \{\; n\}\}$- численное решение, полученное в арифметике конечной точности. Поскольку точное решение $ujn\; \{\backslash \; displaystyle\; u\_\; \{j\}\; ^\; \{n\}\}$должно точно удовлетворять дискретизированному уравнению, ошибка $ϵ\; jn\; \{\backslash \; displaystyle\; \backslash \; epsilon\; \_\; \{j\}\; ^\; \{n\}\}$также должно удовлетворять дискретизированному уравнению. Здесь мы предположили, что $N\; j\; n\; \{\backslash \; displaystyle\; N\_\; \{j\}\; ^\; \{n\}\}$тоже удовлетворяет уравнению (это верно только для машинной точности). Таким образом,

$(2)\; ϵ\; jn\; +\; 1\; =\; ϵ\; jn\; +\; r\; (ϵ\; j\; +\; 1\; n\; -\; 2\; ϵ\; jn\; +\; ϵ\; j\; -\; 1\; n)\; \{\backslash \; displaystyle\; \backslash \; quad\; (2)\; \backslash \; qquad\; \backslash \; epsilon\; \_\; \{j\}\; ^\; \{n\; +\; 1\}\; =\; \backslash \; epsilon\; \_\; \{j\}\; ^\; \{n\}\; +\; r\; \backslash \; left\; (\backslash \; epsilon\; \_\; \{j\; +\; 1\}\; ^\; \{n\}\; -2\; \backslash \; epsilon\; \_\; \{j\}\; ^\; \{n\}\; +\; \backslash \; epsilon\; \_\; \{\; j-1\}\; ^\; \{n\}\; \backslash \; right)\}$

- отношение повторения ошибки. Уравнения (1) и (2) показывают, что как ошибка, так и численное решение имеют одинаковый характер роста или убывания во времени. Для линейных дифференциальных уравнений с периодическим граничным условием, пространственное изменение ошибки может быть расширено в конечный ряд Фурье относительно $x\; \{\backslash \; displaystyle\; x\}$в интервале $L\; \{\backslash \; displaystyle\; L\}$, как

$(3)\; ϵ\; (x,\; t)\; =\; \sum \; m\; =\; -\; MME\; m\; (t)\; eikmx\; \{\backslash \; displaystyle\; \backslash \; quad\; (3)\; \backslash \; qquad\; \backslash \; epsilon\; (x,\; t)\; =\; \backslash \; sum\; \_\; \{m\; =\; -M\}\; ^\; \{M\}\; E\_\; \{m\}\; (t)\; e\; ^\; \{\{i\}\; k\_\; \{m\}\; x\}\}$

где число волны $km\; =\; \pi \; m\; L\; \{\backslash \; displaystyle\; k\_\; \{m\}\; =\; \{\backslash \; frac\; \{\backslash \; pi\; m\}\; \{L\}\}\}$с $m\; =\; -\; M,\dots ,\; -\; 2,\; -\; 1,\; 0,\; 1,\; 2,\dots ,\; M\; \{\backslash \; displaystyle\; m\; =\; -M,\; \backslash \; dots,\; -2,\; -1,0,1,2,\; \backslash \; dots,\; M\}$и $M\; =\; L\; /\; \Delta \; Икс\; \{\backslash \; Displaystyle\; M\; =\; L\; /\; \backslash \; Delta\; x\}$. Зависимость ошибки от времени включается при условии, что амплитуда ошибки $E\; m\; \{\backslash \; displaystyle\; E\_\; \{m\}\}$является функцией времени. Часто делается предположение, что ошибка экспоненциально растет или спадает со временем, но это не обязательно для анализа устойчивости.

Если граничное условие не является периодическим, то мы можем использовать конечный интеграл Фурье относительно $x\; \{\backslash \; displaystyle\; x\}$:

$(4)\; ϵ\; (x,\; t)\; =\; \int \; k\; =\; -\; \pi \; \Delta \; Икс\; знак\; равно\; \pi \; \Delta \; Икс\; Е\; К\; (T)\; eikxdk\; \{\backslash \; Displaystyle\; \backslash \; quad\; (4)\; \backslash \; qquad\; \backslash \; epsilon\; (x,\; t)\; =\; \backslash \; int\; \_\; \{k\; =\; -\; \{\backslash \; frac\; \{\backslash \; pi\}\; \{\backslash \; Delta\; x\; \}\}\}\; ^\; \{k\; =\; \{\backslash \; frac\; \{\backslash \; pi\}\; \{\backslash \; Delta\; x\}\}\}\; E\_\; \{k\}\; (t)\; e\; ^\; \{\{i\}\; kx\}\; dk\}$

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

$(5\; a)\; ϵ\; m\; (x,\; t)\; =\; E\; m\; (t)\; eikmx\; \{\backslash \; displaystyle\; \backslash \; quad\; (5a)\; \backslash \; qquad\; \backslash \; epsilon\; \_\; \{m\}\; (x,\; t)\; =\; E\_\; \{m\}\; (t)\; e\; ^\; \{ik\_\; \{m\}\; x\}\}$

, если используется ряд Фурье или

$(5\; б)\; ϵ\; К\; (Икс,\; T)\; знак\; равно\; Е\; К\; (T)\; eikx\; \{\backslash \; Displaystyle\; \backslash \; quad\; (5b)\; \backslash \; qquad\; \backslash \; epsilon\; \_\; \{k\}\; (x,\; t)\; =\; E\_\; \{k\}\; (t)\; e\; ^\; \{ikx\}\}$

, если используется интеграл Фурье.

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

Характеристики устойчивости могут быть изучены с использованием только этой формы для ошибки без потери общности. Чтобы узнать, как ошибка изменяется с шагом во времени, подставьте уравнение (5b) в уравнение (2), отметив, что

$ϵ\; jn\; =\; E\; m\; (t)\; eikmx\; ϵ\; jn\; +\; 1\; =\; E\; m\; (t\; +\; \Delta \; t)\; eikmx\; ϵ\; J\; +\; 1\; N\; знак\; равно\; E\; m\; (t)\; eikm\; (x\; +\; \Delta \; x)\; ϵ\; j\; -\; 1\; n\; =\; E\; m\; (t)\; eikm\; (x\; -\; \Delta \; x),\; \{\backslash \; displaystyle\; \{\backslash \; begin\; \{align\}\; \backslash \; epsilon\; \_\; \{j\}\; ^\; \{n\}\; =\; E\_\; \{m\}\; (t)\; e\; ^\; \{ik\_\; \{m\}\; x\}\; \backslash \backslash \backslash \; epsilon\; \_\; \{j\}\; ^\; \{n\; +\; 1\}\; =\; E\_\; \{m\}\; (t\; +\; \backslash \; Delta\; t)\; e\; ^\; \{ik\_\; \{m\}\; x\}\; \backslash \backslash \backslash \; epsilon\; \_\; \{j\; +\; 1\}\; ^\; \{n\}\; =\; E\_\; \{m\}\; (t)\; e\; ^\; \{ik\_\; \{m\}\; (x\; +\; \backslash \; Delta\; x)\}\; \backslash \backslash \backslash \; epsilon\; \_\; \{j-1\}\; ^\; \{n\}\; =\; E\_\; \{m\}\; (t)\; e\; ^\; \{ik\_\; \{m\}\; (x-\; \backslash \; Delta\; x)\},\; \backslash \; end\; \{align\}\}\}$

для получения (после упрощение)

$(6)\; E\; m\; (t\; +\; \Delta \; t)\; E\; m\; (t)\; =\; 1\; +\; r\; (eikm\; \Delta \; x\; +\; e\; -\; ikm\; \Delta \; x\; -\; 2).\; \{\backslash \; displaystyle\; \backslash \; quad\; (6)\; \backslash \; qquad\; \{\backslash \; frac\; \{E\_\; \{m\}\; (t\; +\; \backslash \; Delta\; t)\}\; \{E\_\; \{m\}\; (t)\}\}\; =\; 1\; +\; r\; \backslash \; left\; (e\; ^\; \{ik\_\; \{m\}\; \backslash \; Delta\; x\}\; +\; e\; ^\; \{-\; ik\_\; \{m\}\; \backslash \; Delta\; x\}\; -2\; \backslash \; right).\}$

Представляем $\theta \; =\; км\; \Delta \; x\; \in \; [-\; \pi ,\; \pi ]\; \{\backslash \; displaystyle\; \backslash \; theta\; =\; k\_\; \{\; m\}\; \backslash \; Delta\; x\; \backslash \; in\; [-\; \backslash \; pi,\; \backslash \; pi]\}$и используя тождества

$sin\; \; (\theta \; 2)\; =\; ei\; \theta \; /\; 2\; -\; e\; -\; i\; \theta \; /\; 2\; 2\; i\; \to \; грех\; 2\; \; (\theta \; 2)\; знак\; равно\; -\; ei\; \theta \; +\; е\; -\; я\; \theta \; -\; 2\; 4\; \{\backslash \; displaystyle\; \backslash \; qquad\; \backslash \; sin\; \backslash \; left\; (\{\backslash \; frac\; \{\backslash \; theta\}\; \{2\}\}\; \backslash \; right)\; =\; \{\backslash \; frac\; \{e\; ^\; \{i\; \backslash \; theta\; /\; 2\}\; -e\; ^\; \{-\; i\; \backslash \; theta\; /\; 2\}\}\; \{2i\}\}\; \backslash \; qquad\; \backslash \; rightarrow\; \backslash \; qquad\; \backslash \; sin\; ^\; \{2\}\; \backslash \; left\; (\{\backslash \; frac\; \{\backslash \; theta\}\; \{2\}\}\; \backslash \; right)\; =\; -\; \{\backslash \; frac\; \{e\; ^\; \{i\; \backslash \; theta\}\; +\; e\; ^\; \{-\; i\; \backslash \; theta\}\; -2\}\; \{4\}\}\}$

уравнение (6) может быть записано как

$(7)\; E\; м\; (T\; +\; \Delta \; T)\; Е\; м\; (T)\; знак\; равно\; 1–4\; р\; грех\; 2\; \; (\theta \; /\; 2)\; \{\backslash \; Displaystyle\; \backslash \; quad\; (7)\; \backslash \; qquad\; \{\backslash \; frac\; \{E\_\; \{m\}\; (t\; +\; \backslash \; Delta\; t)\}\; \{E\_\; \{m\}\; (t)\}\}\; =\; 1-4r\; \backslash \; sin\; ^\; \{2\}\; (\backslash \; theta\; /\; 2)\}$

Определите коэффициент усиления

$(8)\; G\; \equiv \; E\; m\; (t\; +\; \Delta \; t)\; E\; m\; (t)\; \{\backslash \; displaystyle\; \backslash \; quad\; (8)\; \backslash \; qquad\; G\; \backslash \; Equiv\; \{\backslash \; frac\; \{E\_\; \{m\}\; (t\; +\; \backslash \; Delta\; t)\}\; \{E\_\; \{m\}\; (t)\}\}\}$

Необходимые и достаточное условие, чтобы ошибка оставалась ограниченной это что $|\; G\; |\; \le \; 1.\; \{\backslash \; displaystyle\; \backslash \; vert\; G\; \backslash \; vert\; \backslash \; leq\; 1.\}$Таким образом, из уравнений (7) и (8) условие устойчивости определяется выражением

$(9)\; |\; 1\; -\; 4\; r\; sin\; 2\; \; (\theta \; /\; 2)\; |\; \le \; 1\; \{\backslash \; displaystyle\; \backslash \; quad\; (9)\; \backslash \; qquad\; \backslash \; left\; \backslash \; vert\; 1-4r\; \backslash \; sin\; ^\; \{2\}\; (\backslash \; theta\; /\; 2)\; \backslash \; right\; \backslash \; vert\; \backslash \; leq\; 1\}$

Обратите внимание, что термин $4\; р\; грех\; 2\; \; (\theta \; x\; /\; 2)\; \{\backslash \; displaystyle\; 4r\; \backslash \; sin\; ^\; \{2\}\; (\backslash \; theta\; x\; /\; 2)\}$всегда положительно. Таким образом, чтобы удовлетворить уравнению (9):

$(10)\; 4\; r\; sin\; 2\; \; (\theta \; /\; 2)\; \le \; 2\; \{\backslash \; displaystyle\; \backslash \; quad\; (10)\; \backslash \; qquad\; 4r\; \backslash \; sin\; ^\; \{2\}\; (\backslash \; theta\; /\; 2)\; \backslash \; leq\; 2\}$

Для выполнения вышеуказанного условия для всех $m\; \{\backslash \; displaystyle\; m\}$(и, следовательно, для всех $sin\; 2\; \; (\theta \; x\; /\; 2)\; \{\backslash \; displaystyle\; \backslash \; грех\; ^\; \{2\}\; (\backslash \; theta\; x\; /\; 2)\}$). Наибольшее значение, которое может принимать синусоидальный член, равно 1, и для этого конкретного выбора, если выполняется условие верхнего порога, то так будет и для всех узлов сетки, таким образом, мы имеем

$(11)\; r\; =\; \alpha \; \Delta \; t\; (\Delta \; x)\; 2\; \le \; 1\; 2\; \{\backslash \; displaystyle\; \backslash \; quad\; (11)\; \backslash \; qquad\; r\; =\; \{\backslash \; frac\; \{\backslash \; alpha\; \backslash \; Delta\; t\}\; \{\backslash \; left\; (\backslash \; Delta\; x\; \backslash \; right)\; ^\; \{2\}\}\}\; \backslash \; leq\; \{\backslash \; frac\; \{1\}\; \{2\}\}\}$

Уравнение (11) дает требование устойчивости для схемы FTCS применительно к одномерному уравнению теплопроводности. Он говорит, что для заданного $\Delta \; x\; \{\backslash \; displaystyle\; \backslash \; Delta\; x\}$допустимое значение $\Delta \; t\; \{\backslash \; displaystyle\; \backslash \; Delta\; t\}$должно быть небольшим достаточно, чтобы удовлетворить уравнению (10).

Подобный анализ показывает, что схема FTCS для линейной адвекции безусловно нестабильна.

Ссылки